{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extended Exercise\n",
    "\n",
    "Download the Quandl dataset \"Capital Market Yields – Government Bonds - Daily\" from Quandl (code `RBA/F02`), and perform the following steps:\n",
    "\n",
    "1. Split the data into a training set and a testing set. The test set contains the last row of the data. The training set contains all the previous data.\n",
    "2. *With only the training dataset*, fit a PCA model and extract the principle components\n",
    "3. How many principle components are needed to fit at least 99% of the variance in the data?\n",
    "4. Fit a OLS estimator on the lagged data, *for each principle component*. I explain this further below.\n",
    "5. Project the data into the test period, by taking the bias and coefficients of the OLS and predicting the next values.\n",
    "6. Map the projections onto the original data space, by taking the PCA's inverse, to get the actual predicted values\n",
    "7. Calculate the sum of squared error for your predictions.\n",
    "\n",
    "\n",
    "\n",
    "In step 4 above, I ask that we fit an OLS on each principle component. If we have a principle component $C$, we want to use the value $c_{n-1}$ to predict the value $c_{n}$ using the linear equation $c_n = \\beta_0 + \\beta_1 c_{n-1}$. In other words, we fit a linear model saying \"given the a value of this variable in a time series, what is the next value?\". This type of analysis, moving average, will be reviewed further in a future notebook (1.6.1). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: This exercise fails for some versions of statsmodels, including different computers running 0.9.0. There is a simple fix, noted below where the error is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR:root:File `'setup.ipy.py'` not found.\n"
     ]
    }
   ],
   "source": [
    "%run setup.ipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "import quandl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bonds = quandl.get(\"RBA/F02\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bonds = bonds[bonds.columns[:4]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bonds.columns = [\"2Year\", \"3Year\", \"5Year\", \"10Year\"]  # I plan to use statsmodels, so no spaces in variable names!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>2Year</th>\n",
       "      <th>3Year</th>\n",
       "      <th>5Year</th>\n",
       "      <th>10Year</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2013-05-20</th>\n",
       "      <td>2.480</td>\n",
       "      <td>2.530</td>\n",
       "      <td>2.760</td>\n",
       "      <td>3.220</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-05-21</th>\n",
       "      <td>2.520</td>\n",
       "      <td>2.570</td>\n",
       "      <td>2.795</td>\n",
       "      <td>3.255</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-05-22</th>\n",
       "      <td>2.535</td>\n",
       "      <td>2.580</td>\n",
       "      <td>2.805</td>\n",
       "      <td>3.265</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-05-23</th>\n",
       "      <td>2.530</td>\n",
       "      <td>2.585</td>\n",
       "      <td>2.825</td>\n",
       "      <td>3.295</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2013-05-24</th>\n",
       "      <td>2.540</td>\n",
       "      <td>2.595</td>\n",
       "      <td>2.840</td>\n",
       "      <td>3.310</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            2Year  3Year  5Year  10Year\n",
       "Date                                   \n",
       "2013-05-20  2.480  2.530  2.760   3.220\n",
       "2013-05-21  2.520  2.570  2.795   3.255\n",
       "2013-05-22  2.535  2.580  2.805   3.265\n",
       "2013-05-23  2.530  2.585  2.825   3.295\n",
       "2013-05-24  2.540  2.595  2.840   3.310"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bonds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f5efb78d9b0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bonds.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "training = bonds.iloc[:len(bonds)-1].copy()\n",
    "testing = bonds.iloc[len(bonds)-1:].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1514, 1, 1515)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(training), len(testing), len(bonds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DatetimeIndex([], dtype='datetime64[ns]', name='Date', freq=None)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Confirm we have no overlap\n",
    "training.index.intersection(testing.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Confirm we didn't miss anything\n",
    "(training.index.union(testing.index) == bonds.index).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PCA()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,\n",
       "  svd_solver='auto', tol=0.0, whiten=False)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.fit(training)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([9.84317967e-01, 1.25451480e-02, 2.98756705e-03, 1.49317754e-04])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.explained_variance_ratio_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Note - only two features are needed!\n",
    "9.8e-01 > 0.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(copy=True, iterated_power='auto', n_components=0.99, random_state=None,\n",
       "  svd_solver='auto', tol=0.0, whiten=False)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# You can simplify the above using this parameter:\n",
    "model = PCA(n_components=0.99)  # i.e. explain 95% or more\n",
    "model.fit(training)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_training_proj = model.transform(training)  # Project the data onto the components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1514, 2)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_training_proj.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5ef80feac8>,\n",
       " <matplotlib.lines.Line2D at 0x7f5ef80fec18>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXd4FFX3x793e3onhJrQO9JBQJoUu6JSbFix62svvPb66qvYC/afL6Io2EBFUQRBpPfeIRAgIaSXbfP7487szOzObjbJ9pzP8+TZmTtl705mzpx77ilMEAQQBEEQsYMu3B0gCIIgAgsJdoIgiBiDBDtBEESMQYKdIAgixiDBThAEEWOQYCcIgogxSLATBEHEGCTYCYIgYgwS7ARBEDGGIRxfmpmZKeTm5objqwmCIKKWdevWFQmCkFXXfmER7Lm5uVi7dm04vpogCCJqYYwd8mc/MsUQBEHEGCTYCYIgYgwS7ARBEDEGCXaCIIgYgwQ7QRBEjEGCnSAIIsYgwU4QBBFjkGCvg1+2HsfJ8ppwd4MgCMJvSLD7oMbmwC3/W4drPlod7q4QBEH4DQl2H5TV2AAAR4qrwtwTgiAI/yHB7gVBEFBeYwcA6HUMAPDTlgLkPrzQJfAJgiAikbDkiolkft5SgO82HsWibSfQOj0OAFAmCvjXFu8GALz9x148cm7XsPWRIAjCF6Sxu3Hr7PVYtO0EAOBIcbVqW5XVAQB4f9l+v8+3dHchDp2qDFwHCYIg6qDJC3ZBEOB0Cn7tm3+6WrPd7nB6bZ/28Wpc9PaKBvePIAiivjR5wT7zt91o9+hPsNq1hbNErd2h2b6vsAIdZvyMeevyUWOT9xEEAcdKuJtkSRXZ5AmCCB1NSrBvOlKC5xZuhyDIGvobf+wFANTYHdhzotzrsZW1Dpj0npdrc34JAOC+rzehz9O/udqf+nE7znp5SaC6ThAE4TeNFuyMsdaMsSWMse2MsW2MsbsD0bFgcNHbK/DBXwdgc3iaXuwOATd85ln8wyB6xGzKL4FVYXKRNHybXT5Xtc2BwvJavLxoJz79+2CAe08QBOEfgdDY7QDuEwShG4DBAG5njHULwHmDhkPDpm5zODXb7WLb64v3qNr/++suAECtm3398vf+xttL9gWqqwRBEPWm0e6OgiAUACgQl8sZYzsAtASwvbHnDiSHT8lBRjanE3HQq7Zb7U7ofLzmNh4pUa3PWrYf7bMSUFVrV7UfPOUZzCT5wRMEQYSCgPqxM8ZyAfQBsEpj23QA0wGgTZs2gfxavyitlicwHaIp5u99Ra42u1PwcG+si4fmbfFrP4dTQJXVjngThQ0QBBF8AjZ5yhhLBDAPwL8EQShz3y4IwixBEPoLgtA/K6vOItsBx6mYMLU5uflk/aHTcpsXl8WrB7dVrbdMjWvQ95/3xvIGHUcQBFFfAiLYGWNGcKE+WxCE+YE4p7+UVtv8CvGvVbgz2kWNPSXO6Gpzd3dMsnDt+sbhea62KQNaY9mDo1T7tUrzT9AfKKIgJYIgQkMgvGIYgI8A7BAE4dXGd6l+9H7qV/R68tc696u0yrbwB77ZBABIMMumEbti4jQz0YT/XNoLZoMOzZIsrvbrhuZBr2PISjK72rQmXJW8PuUMAMBZnUI/SiEIomkSCI19KICrAYxmjG0U/84NwHkDSlWtHDy0Yu8pOJyCS3MHgANFFWDiHOdbV/TFuT1zsOvZcxBnkidZ48Vl5VxotU07cAkABuam46IzWmJgXjpW7ivCir1FXvclCIIIFI0W7IIgLBcEgQmC0EsQhDPEv58C0blAotTYAcDudKqiSe/5ahMEAbhvbCcMbpeh2veZi3ugX9s0l31dx2TJXm114KxOWfjkugEe3/n8xJ4AgCSzATaHgCs/9JhTJgiCCDhNxk3D3S3R4RRUdncJi1Hv0Xb14LaqSVSlYK+1O9G3TSp6tkzxOK59VgIAYMmukw3uN0EQRH1pMikFKq1qk4ndKeCN3/d47GcxeQp2d9z90uOMepgMnpeSiS8AP3OMEQRBBIQmI9ir3QS70ym48qxnJJhc7RYNAe2Oh2A36TXzyEj46zlDEAQRCGJesAuCgCPFVai2OZDoxQvmkj4tXctaphh3mFsgaZxRLdh7tkzBw+d0ca3Pu/VM1/LKfafq1X+CIIj6EvOC/dsNRzH8pSX4c9dJldB2OgVYjDrcNDxPZUaJ80Owu/u8x5n00Cm0+FtHtsctI9q71rOTLXj+Ej6ROvWDfxr8WwiCIPwh5gX7LjEV777CSsSb9HhR9FTZe7ICNTYnLG72cX809hqbWrBnJppV68qRgURLhTmm0m0ilyAIIpDEvGDPTJCFbpxR77KPXyG6HnoK9rovyVmdMlXrOSkW1XqChmDPUgj/YyX1y0lDEARRH2JesOeflrMtWkx6GPRqA7nZoMPA3HR5Hz809hcn9lKtZyerBbuWxt41JwkX9m4h9okEO0EQwSPmBftnKw+5luOMOjjdXNctRj36tU1TrdeFyaBTe9K4HZNg9jwHYwwzzusKAMgnjZ0giCAS84JdSbzJgCq3FABxRj0YYzCKmrw/phgAuFIMWDLqPXOta2nsADfHmPQ6rDlQjHJF4rJftx3H1qOlfn0vQRBEXcRM5Omy3YUQAIxQJNtyukUGxRn1HhGokrZtMephc9hhNtStsQPAPWd3RJ/Wqeifm+axTcvGDgA6HYPV4cQPm47hQFElfrxzGABg+ufrAAAHXzzPr+8mCILwRcxo7Nd8vBrTPl6taquxq7Vzi1HvEYEqaeg3DOPpeb1p2+4wxjCqSzMkWeTUv9K5jD6ClSS2kIZOEESQiBmNXQv3aNM4k87DB13S2O8e0xG3j+rgl1D2xqJ/nYUdBeUNPp4gCCIQxIzGrkWVm2CPNxlwy4h2qrbsZO6GyO3sjbscbTMSMKFHc5/7SFGoUqZIb5Wb/GHtwWJc8OZy1PhIHUwQRNMjpgW7e650i1GP1HiTKlNjh2ZJIe1Tv7ZpuGl4HoorrQCAP3Y2PPPjfV9vwpajpThc7FlAmyCIpktMC3Z3jX2baNdW5nEJB5mJZlTbHJj8/krcLE6cAp6TvXUhvRx0no45BEE0YaJasGsJQmWpOsnGLplbpNqo/uSDCSYZYhTqqgPFqvYKa/1SDUi/3z3FAUEQTZuoFux2DcFepRCO1Ta+/PYVfXFh7xZ49mKeJ0YXZhU3I9Gk2V5aVXdRbiXSz7c2wk5PEETsEdVeMU5BS7A7XC6I1VYu8JIsRrwxtU9I++YLZf4aJaXVNrSux3mk3+/u6UMQRNMm5jT2CkUAkuQtomV6ef6SnvjipkHB65wP0hKMqvVmSVzQ3zVng9/nEAS5tJ9WiT+CIJouUa2xOzQEe63C3iwFKGmlCbhiUJvgdawOWqXFq9ZbpMbhZHkt9hdV+n2OX7efcC2Txk4QhJKo1ti1Jk+V9mZp8tQc5snSupB82uvD+0v3uZZr7eTHThCETFQLdi1TzHMLt8PuUJso/E3sFS66tUh2LZ8sr/HrmJQ42Zwzf/3RgPeJIIjoJbIlXh1oTZ6uOXgaD83bglMVtaixOaBj8FloOhI4s32Ga/mHjcfqfXxjgpwIgog9Ilvi1YGksU8dqLaXz1ufj+s/XYMamwMWMS1vpDEwjxf3OLdnc/RqlepqT44zejtERVmNHV2ahzZqliCI6CCqJ08lG3u/tmm49sxcjH9tmWvbgaJK9GqV6lfhjHAw9+Yhmu0GP3zsBUHAjoIynN01GzuPU9IxgiDURLXGLnnF6HVqmzPAg3eqbQ5YDNHxE2ffyF0v/Sl0bXcKqLI6kJuZEOxuEQQRhUSH1POC3SXYdWjuVlC6otaOiho7LKbI1Njd6dOGm2Pc88VrIWWEjI+S30YQRGiJasEuTZ7qRRv685f0VG3/ZdtxWPysiBRueIk+eFR40kLyWzfpdRjeMRMAXJ5ABEEQUS3Y7Q7ZFAMAUwd6BuRHuqujBGMMcUY93vhjL3YeL1Ntq7E5VHnbJV99k0GHYR24YKfoU4IgJKJD6nnBpbHr+M/Q8n7xVn80EpHSDF/3yRpX2+oDxejy2C8Y++pSV5tNfKGZ9DrX5DAJdoIgJKJasCsnT70R7hS9DaGgtAYr9hYBAJ77aQcA4OApuZiGzS5r7Mlx/MW15wR5xxAEwYlqwa6cPJUwuXnBRNME4xmtZX/22asOAQD0ikFItdWB5XuK8N1GHmlaVmNDvzbcH35vYUXoOkoQREQT1YLdffIUAN67qq9qn0RL9JhilLVLpfwxhRW1rrb7v96Eqz5ahXf+5HliiiutyEnl3kBzVh8OYU8Jgohkolqwy5OnsmAf3SUbE/u2dK13bp7scVykct3QXNdytc2BKqsdR4qrAQCds5Pwz/5TAGSvmFtGtHcV4N56tIyKWhMEASAKBbsgCPi/lQdRbXUoJk/Vk6bn9cxxLbdMVfu3RzKTB8ipESprHSirll0fE8x6GPTq3+keVeteao8giKZJ1An2JbtO4vHvt+GFn3d4nTxVTpgmmv3LvRIp3Dm6AwAeYPX491td7XtOVOBEmWyWGZib7nGsPz7wBEHEPgER7IyxjxljJxljW+veu3FIZogfNh1TCHb1zxjUTs6WmJMSPRo7ANw3rjP6t01DldWuKqZR7ia0O2QnehxbXGUNev8Igoh8AqWxfwpgQoDO5RPJ66WkyobtBTyQR+/mv67XMex//lysmXE2WqfHe5wj0kkwG1BR60B3MU97qzTPQhzKUYkUmHWi1L9c7gRBxDYBEeyCICwDEBIDr1Fhd3nl110APG3sAKDTMWQlaReNjnQSzHpU1tpRUmXDxL4t0UVjAljpxvnCxF7ITjajgAQ7QRCIQhu7smiStKwl2KOZOKMBe09WoLzGhkSzAac1TCzu2Syzky04UV7rsR9BEE2PkAl2xth0xthaxtjawsLCBp9Hq3BzhBdIqjdLdvGKSGU1dhj1OpRV2zz2SYs3qdZT4oya+xEE0fQImUgUBGGWIAj9BUHon5WV1eDz2DSyGOoisEJSY1AOQFbsLYJnAUAgPUEt2JMsBmw8UoJdisIbxZVWbMkvDVIvCYKIVKJO19US7O5miWjHrEg1/MD4zpr7ZCS6CXbRrXP8a8uw+kAxxr66FH2f+Q0XvLVcc5RDEETsEih3xzkAVgLozBjLZ4zdEIjzalFZ6xld6W6WiHauH5bnWu7YLAmComj3qM5ZuLB3C3TNUU+oZivcOie9vxJ7Tsq5Yw4XV4EgiKZDQBKpCIIwNRDn8YdHv93i0aaLscnTG4blYeHmY1h/uARxCu+XYR0y8cbUPkiyeI5QOmn4tUuU15DtnSCaEtGTIQtQaa4S156ZG/qOhIBZ1/THX3sKkZVkdtnYH7+gm6ZQB4DMRO+unSVVJNgJoikRVTZ2ZUi9xBMXdAtDT4JPZqIZl/RppWrzNUnczkdh65v/ty5g/SIIIvKJKsG+S6OYhFbVpKZIs2TvqRNo8pQgmhZRJdjzT/NJwCdjVEv3ipa/o59cdEYLxBn1qKAEYQTRZIgqwV5ew4XT5f09i1bHMmmiz7pRX//RyWX9WqHa5sDGwyWB7hZBEBFKVE2eltfYoNexqCp3FwjevbIvftpSgLYZ3u3oWtw8oh3aZXFvmSOnyeWRIJoKUSbY7UiyGJqcXb1ZsgXXDs2re0cFW58aj0SzAWWiq2MlmWIIoskQVYL9ojNaok+b1Lp3JJBo5v/aRBP/LKshwU4QTYWosrH3a5vm4QJIyLw+5QyPNp2OIdFsQAUJdoJoMkSVxq5k2pC2MBublq29Li46oyVOV1oxIE9dNs9i1KPGToWuCaKpELWC/amLeoS7CxGJli3ebNChxkaCnSCaClFliiEahsWoQy0FKRFEk4EEexPAbNCj1kaCnSCaCiTYmwBmow61ZGMniCYDCfYmgMWgJ3dHgmhCkGBvAnRunoRNR0pwhApuEESTgAR7EyBTLKM3/KUlYe4JQRChgAR7EyDOFLVerQRBNAAS7E2AppY0jSCaOiTYmwAWI/2bCaIpQU98E8Cgo38zQTQl6IlvAigLdBRXWsPYE4IgQgEJ9iaAUmN/8JvNYewJQRChgAR7E8Cg0NhPlNWEsScEQYQCEuxNgFZpca5lyvJIELEPCfYmQIdmSeiakwwAqCbBThAxDwn2JsLH1/YHAAzvmBXmnhAEEWxIsDcRclLikJVkRq3NgdmrDsHpFMLdJYIgggTFmjch4ox6zN9wFPM3HEVGggkTeuSEu0sEQQQB0tj94cQ2YGZPoPRouHvSKJSpBSiNL0HELiTY/eHdM4HSw8BH4/h6TRlw+lB4+9QALMri36GwxNitUXmdCCLaiV3BvnEO8MFoQAigBCvLB17tBvy3I/B6Ly7go4iDpypdy0dLqoPyHa/+thsfL94IFO0Ffn+KX6eV7wTlu2IKJ5UuJAJHbAr2sgLgu1uAo+uAksONP58xQXHuo4BdDPJ5sTVgr238+UNESZXNtfz673sCcs5V+0/htCJNwRu/78EZy24A3uoHbJ3HGxc9Aiz7r+fBexcDc6cB278PSF+iktKjwOcTgafTgI1feG53enFP3fIN8GQK8EpXoGAT8PdbwM6fgLUfA0fWAG/2AzbMBk7tA45tCO5vICKO2Jw8/UshRLZ9CzTrCnQa37BzndoH2CqB4fcD8enA8plAZaG8vWAT0HogX3bYAGsFEJfW8L5HEU6ngMmz/kHXnGT8fPdwV3tf3V6+UF4g7/zHM8DB5cBV8wCdaBL636X8c/t3wAP7gITMEPU8AijaCyQ2A+ZezRUQAPjrVeCMK+R91n4MLLgH6DgeuOht4JeHAIcV2PGjvE/5MeD9s7S/4/vb5OW7NgDp7QL/O4iIJDY1dsgh9Fj8BPDFJODtQUBNaf1Plb+Wf/a8DBhyO9B7inr7PwozwzOZwH9yA2v+iWAc4u/cUSCbpN41zvR+wP4lwP8mam/b+7t2+387Aa/1AhxBnuz9/WmuOYfifycIfETzYmtZqANAXKp6vwX38M89i4D/duAjIKVQ90V6e/X6m/0a3t8QU2t34I4v1mN/YQUAYM+Jcmw92oBntwkTm4K9qgjI6AAMuUNuK9wJ/P1m/c7jsAHfTufLGR3557hngYePyPsUbOKfdkXWxCOrvQ+hI4jG+rI73I9f+TbO0a/hy+1H88/rfgGmLZD32f8nUFkEVBXz9bHP8BHOoeXqc9lqgG9vASpOACWHgD9faFRf6+SvV4B9v/PPYPPzg55tPScBp8SRzqG/vZsQMzoAeWcBN/8FPFkKPHgA6HEp0HEc0G4k8O9CPvq5eSkweTbwuHidBSdQdszzfJVFXOgvvD8QvywgrD9UggWbCzD6laUAgLEzl+H8N5fXcRShJDZNMaX5QGI2MP45Loi3zQe+uZ7b3uvD4X/4Z49LAb3iUpmT5OWSw8AbfYDi/XLbx+OAyz7mx0UQLVPjVJOmVocTFl3DqyvZRcGehdPAwRXAokfljVd8jWfnLkXlugS8MLEX8PhpYN4N/H+x7VugRR++X0YHAAxY/3+AIQ4Yfi83db03TP1l274FxjwGABAEAeW1diRbjA3uu1f+eAbodhGQ2THw5wa4IF09S912/x5gy9fAlrncbq5k/AtARnugugToPdnzfPHp/F5TYhBNWl3P559D7gBWvgUsfgqY+D6wahbw8wPqY07tBfpcKf9f6qLkMJDaxr99lVSXcFOc8hlyw2Rgmu1OpwCrwwmnICCeyj36JHY09uID/KFY+hKQv0a+QRnjAjZ3OFC0y//zHfiLa04AMPrf6m2MAY/kA2OfBpx2tVCXkDT5COKbW4eo1mttjfPEcDgEtGPHsMZyO/Dpua72S2qfwtFyGz7cVIM5q8XRjU4HjHiIL/90P48NALhgl7T71e8Dr3RWC/XMTlwbLd7H/8cAHpm/Bb2e/BXV1gCNigp3q9ff6h+8mAXpXul8HvDIUf6X2Axoe6b2/r2n8PkhLaHuL+OfA0xJwKk9QG2Fp1CXmDVSNj2643Rwc1XxAWDHAuC1nnxEVbi7TgeCfYUVOFpSjTOfXwz8py3wQis+d+UFb/flifIajHz5T3R7fJHP74s4KouAfX/wvxAREMHOGJvAGNvFGNvLGHs4EOesN6ve559LnuOfOWeotye3AE7uBH56EDh9kHsQeKPsGPDZ+cCfz4vHtvTcx5wEdDpHXk/LBe7dIa+veL2+vyDo5KTEqdYLKxrn0bPou//DH2bPIfwGoSOGvqhxEyu14B/v4p9pbYGJs4BuF3vuf+1C4I41QNcLAABvvvokch9eiC/X8JdFeY2Na4DzbgSK9nD3U2tV/X/I5i/559lPyW3zbtA2XTQGWzWwcyFfHvEAYE7kfwBXRIbeDfSawk0sfacBl8ziGnkg6HIuf6ms+ZCvJzaXt804DowUR1sfjvG041srgafTuZnqi0l8ZAEAm+YAbw8Anm2mNkW6MeaVpRj64h8YVvmL3PhmXz5y0KCiVp5PERRzHgeLqnA82tJO22uBl9sDn1/C/zZ/Hfz5IgRAsDPG9ADeBnAOgG4ApjLGujX2vPXGqBZaHl4w8RmAtZxrha/3Br6cCvyk0FyqTwO7f+UTW692VR9rMGt/Z1Yn4IkS4F9bgbs38ZfHjYpJwI1zIm4i9cc7hmHKgNYAgPGvLfO+Y02Zz4cVRXsxafe9Hs3TrA95P0anB25aIq+f/RS/tjo9MOkzYMBNfGQ14zi/rrmi5t7nagDAnYbvEAf5wS6rscO5exE3Y3x5BfDOEOD5HD5yezIF2CUKkhPb+AvAG7UV3Aw07F/ADb8Bxnjg8ErP+0ALQZDnC+riuebAitf4cmpbz+1jn+amEgC48A2vWrrd4YTVXs/RVpfz+T2++AmgZT/gvp1As+5AXDp/dkY+BNz2D//tX13Fr1/lKX6s0sRWtJu7p7YaqD7/s1nArFFeX6wMTrxk/EDduPgJzX2/XpfvWi6rloXgIUUcRtTg7rAx/0Zgxw9B/9pAaOwDAewVBGG/IAhWAF8CuCgA5/WfFW8Ay1+V17tdDFiS1ftkuHkJANzWWVUM/PY48O5Q4IvL1V4uV36jFkRaMAaktpbXW/UHLv+UL393C/f2KM3XPDQc9GyVgjM7cBusx+QnwAXgTw9wj41ns7g3kd0KHN8KrP5ADqTZ9q3qsI3nL8RY6ytY6uztuwMt+3KN9MlSLkiVnPdf4NoFXNAwhZ1VsbzDcj3G6Lgnyd0zP4NOmtwu2s0DyJTMmcxNC++eCfxwB7xiq5RdVFsPBG792/dvAGTT31OpwEt5wM8+XmiAWkvrNaVRmvil761Ep3//XL+Dul4ADJzOR5YDp/Nreu0C4A6F6aVZV+AaRUzBe8P4b1z3KV+f9H/ytuzufAJ3vGJS+9h6eSTmRk92wLV8n/UWvmCrAt45EyiUTaSCIOC37Sdc65/+fdC1fEAh2I+XRonmXluuXs/uAXS/JOhfGwjB3hKAwk0E+WJbaHDYgd8eU7f10tB0+l0vL5/7X3li86U8bjYpE22qix7lk0KPnQI6juWCqL50uxgYeDNfnn0pMLM7sHsRsPJt31pwiDDpffzbf39aPblXuJML+PeGctv4s1l8CO2UBdXltY/j4m9KsccZvKRio2vl2ISPTK+gOzuAheZH1TvpTcA5L/HJxOwevO3DMfyz6rT3k5ccAQwmeT09T374Nn/Ng3+q3Y531zZXvcdfdk4nn6RUCCuUHQOeyRB/yGOyVt5ANh3xMfrwBmPAuS/zkaXkshufDiRkqPdrPRC4Zzs3DZXLpqiD4z9FbafzgdtXA4NuAcY+BeT0AobcBjx0ELj8M77jlq+BE9s9vv4Bw1cAgOnWezDPeRYm14rP7Mlt3ONJHNm65zCauXg3dHAiEVU4VCSPBhZuqacjRLg4KV6L817hcufmv9RKS5AI2eQpY2w6Y2wtY2xtYWFh3Qf4y+mD6vXHTnF7ojs6HXe7mzwbGHgTFwDMy88fcofaC6a+MAac+xLXjiS+mMRfGl6Gny6czqALf4vRy+8+sgZY+5Hvg512/jsqjqNISEZuzWysEbr4POSf/acaXZJvv9ACT9quca0vNM8AACxz9ORmmydLgccKgUE385f21d+pT2DzYnuvOAkcWOp5H3W9kH/OvxH4dQbwxRTP4yRmHAfS8oCvr+XXZvmrwNsKU8XxLfJyel7dPzbcpLQE+orX+sw7sWzyNoz83oSnf9wOZHUGzvkPYFF478SlAd0vBgaLAVFbv3Ftql79fzhouQLD9VsBAL85uT/9KqEL8ke/xXfa9i2whM9nSZr4tCGyqWqW8RVstdwIw0n5Ojoak4Lh1D7/4wEayvGtfLTz1VV8PecMLnd0oRG5gfiWowAUtgi0EttUCIIwSxCE/oIg9M/KCmCxh6oieXnko74Fct5w2QUsIROY/qf2floaf0O4ayN381PyzzvA8tc89y3YDHx4Ns9D83wL2b4ZBOKMXlwcC3eq1x/xYkKKTwdObEO+kAllMJi3kcCUWf9g2serG9BTNZ86JqBnzYeqtift07Q1oMQsYNLnwMhHuMA5tl62hQsCn0isLQcOreBtg25VH+8+mXvkH35cRSEwZyq3wQPAPdu46UiKGF31rnzM8a3AoZV8klEiPnDRtaWKFBEBp+804LpfMMtyHa75jHt4LdtTh0I24QWgeS856KriJOJ+utO1+UnbNRBcIodhdeIonFsrOiiIdueCUu6O27NVKoyw42HDHJyt5ykR3iq/G0Zwjb6mMR5d39/OBa4ilYXN4URZTQCup9PJrQhfKiKI80bweY0QEghn0DUAOjLG8sAF+hQAV/g+JIBI4f03Lam/2aR5L+DCt7g7nSmBRzmmtPKMAGwojPG/6X/yhGTJrXiWyMVPcE+DnQv40NxpA+Zeoz725XbcVh8Ee1zv1vz3GfVuAlHKgXP/Xq6F6Q18nmHpf7h9VW/m/ao4AVScQC5LUB1udXh/2A4XN8BbRaRYkYumHPHoUfMhtlpu5Nv0PgRltwv538Hl/IX6koa2fMEb/HPI7ep2nQ64+luuyUsRoH88wyM6d4keVaNm8PsFAIbdI3tkSbw3VHE+AzD8PnlCOAD8Z9FOPH90WpK4AAAgAElEQVRJz4CdT4VOD7QdguffXehqOlHmhxdV856yW99/ZS+oZ/uvxKfLD6h2ffCbzbALuXjdPhF3F34L2Kpd+Yz6tknFL9nvoH0pjyU5qW+OZo7jGKjbgRXOniitrr8Qtjmc+N/i1bhOeinPvQaY+AGQ1BwP/noafxUmYPWjY6DTNcBUsuNH4OQOz3ug7TBg6pyQmF+UNFpjFwTBDuAOAIsA7AAwVxCEbY09b52cPgjsWSwPdRqSZ4QxoO/VQFI2dzu7fzdwsw9PkYbSog/wxGngni3AUHHCcNlL3P725VRPoS4xf3rg+wKevndIuwzYHAI+/Evhg192FNAZZaEO8HmGGxdzj5+EDK4BixwQPG3q0kvDHbOh4bfafXM3AgDGdcsGAFQg3rUtJ8uP/7svbUmaLDcleG5rPxrof718T+z4kb+MJZS5V/RGbptu3lN7JNiyHzDqUTlPjgYFpdX18s0PSoCWAvfC5+nxprpHCQlZQHkBnBtmu5rG176Ithny/2x4R/4/kwLcdjpbAxBQeXQbrNXl0MGJ5n/NcAn1G6z34Zu+nwMAZptewFemp2E5vZuPwJa+DFgrcex0JSa8tgzHS2u0nQIAzF+fj6JloqlR8kqafxPw2QWYWXANJlQvwNfrjmge68HGOcCz2cDR9dyM9NVVnkK95+XAdQu1760gE5DwLUEQfgLgwzE8CLzu5n2R1KLx55R8ioPJ2Ke4kNQKK797Ex8eZnbmHjXVJdwVLwj9SrLwf/2zC3dgysA2SDTqeOBJVmff5qwz73SF919t9QxZMLmPAkSMviZs6+CUqLErC3G/ZJuMQqTA5M8Lw90V1pQInHU/sPhJ2fxk8nGNc3oDPS7jtuOi3UDvK/joUPSvd5HSErjFLfT9ht/46KxZ3a6TQ174A20z4rH0gVG+f46eweYQvJvUAkRhuVpDP15Wg/7P/YYJPXJw81nt0KNliudBCdzMqhMTkH2Z9ywKDrTDpf1aYXN+KR6c0AVZSWbkPiyPBHYLfNST8OkYTALQ2tgN8Zv5pOP7vb7GjOFnIv90NbcNABik24lB+6YBXw7hJrElzwKJPXB/qRHNZyoyWd66kucnsqQAm+dCv1fAVYat2BfXE4eHvo0Ru5+Dbo8c7PSs8RO8930RkDCJR533nsIniJU4ncCaD+Tn9wMf/yv3wMYQEp1xue4TJ90ubtxkZ6jJG8E/9WZg8C3AmCflSZWhd/PPiR9yj5qja7kJZ+Vb3MWsx6UNc5Ur3A0s+BefNG7eA90cO7EccaiCBQ/N24y3m/3IIxMH3OT7PKYE4LqfAXMSKl7j+UzGdGmGNQeLUVZjh0ExOTS5f2t8tZZrQCfLa1FcaUV6gknztL7ITORxBEUVsklmU971MOh0KC+uwrMLtuOuszv61mDv38tHRxNe5C6pABfsEoY6+tXhbHlSsP/1QOsBvvfvfz2PXG41wK9huJS359Ap3yYrq90Jm4Pva3UENx9RlcboweYQ8OOmYzhSXIXvbh/qeVD/6/hks8gaZxc0TzEh3mTAy5dru8IeFJqr1ofoRU+Se3fg5mSusJXV2HGV9RH8z6Rwr5RMKgBaVGxFC/f33LvqSOvLxO1Xld6M5fPyAUzDwUdfx69zZuKdgy3xnflx3GL4EZgrTqz+8zafJ6s+zV2i9SbRzOQlNqXbRTyPT59reIoGpfNEiIkiaahAmozqOB649ENPn/VIp1kXbj/PG+FdSLceAIBxH+Lt3/MkTgDP8Hf9L9rH+GLtR3yiULT7/gtAe+Ng3Gm7Cws3F+Ct9M/4NGhfL2YhJWL4u1F/BL1bpeKjawdg/MxlKKsph1GhQceZ1E/aW3/sxeMXaMeuHS+tgdXuRBvFkF2iZSrXuEd2zoLTKSDRYsAn1w7EvXM3Yn9RJfYvPwDGgBnn+YiLS8ziJiUll3/KPVkurcMTCODa28bZPHq2LqEOAOe9yidb/bStVlj9i0asVERlSgI+WFSKfWqdHocjxerCLMlxXl6ipgR8f/F2LJv7BmzQ44ddVozopKHZK3BA7ym0AT6ylRYtBix39sB11gfQoUMnzDgsKiDT/+RmkdXchXSVsws2O9vhhkHZ0K37RPP7ljt7uJbvWFCA9u1uxMYDe/CAbTpu038PpLZGXpno3z/7ckBwaKcNUTL1K6DzBHk9O/QxmkqiT7Cf3Cnnme5zVfQJdYm6JkUtKQAEj0AgHF7JRyySZvzPe9y7o91IHmhjTuIjgR/u4C5WF4geOKve8/iKjow7L/Vi+8CqRC8c96GnD3SMoV9bHtgTb+ZC3KiYeHIX7C1SLV7PNfgFHrF78MXzPLZJ57l/XGc8NEF2rdQrvssfIed0CuqJse6X+D85zRgw7Uf/J8GkiXM/WbXfv+hVZbj98j1FPvZsPFW1XGOfOekMTHp/JZSm6wSTdzNQabUN85xyjnjpHvHFcqc8CfyO/ULcdk5/1faUOCMAhiXOPmibmYsJe16EHk4sbNEHaNEHf6/8CxZmxWTr4wCAIZ07oUfZMT4ibtkfDmMCnnr8HvF75P/Lgs0FuGNUB+h1DMfbXYZRe0YivdKE9U+O5e6KpxQFacY+w2NcEpvxSfPkViFzX6wv0SXY9/wGzL5MXm83Inx9CQXjX+DVhwA+Ouk1iecwWf8Z0GMi8KIiu96SZz2PP7aBaxrjnuVJoKzqKLguuiNIQhWuN9QzilHE5nC6bOdSoWylLT3RrL69Gpq0y+ZwIslsUAlyAMjLlCel6pKhR4qrMPylJXhjah9c2LuB8zFB9GyQJoh9cehUJcoVATzbC4JbmlF6icSbDHjigu544gfZJ6LcRzF09wnWDs085y+uHtwWn/9zCDkpFhSIvutjal8GAPz0zA2AQf3iSFKY2Yx6hp0Cv/cFQUBFrR1X2h6FTmEimTJ7P7Y+NVfub5UV/+fQLraz60Q5kiwGXNKnJf5SvixH/Vt+rh49FpZJ0IYSma8bb+z5VV6+d4c6SCIWGXIbD765aj5w2Udy/psF/1ILdV8cWAq8P1wt1Kd+yXOjAPje9G9crBdD6OtKn6DA6RTgFACDOFkqpVE1KCZPlYLdYtShotaOUxW1WHPQz9wqInaHoDqvxJgu2a5lBt9Cd9dx/vvnr4+c9A5KJOEnTWq7U1heixEv/4kZ327R3B4MqkRTTIJZj2ln5mLH07KpYfneIpRUaQfSlbi5ImYleeZaeubiHtj+9Hi8faXsorxPaIlRQ4fBbPAcDSgnyZVetZvyS/HI/C0QoIMD8nHZyervlNwoc1I8R41bj5YiLd6EiX1bYdqQtiiutHLPmhEP8IDHB/ZFlVAHok2wdxGH6WOeUNnfYhrGgA5juInFnOSZL/u+XfzGU3L3Zv5CuNbNUUkyO2R3d7nwtdMdd222JmskpvKCTZzAljR0aWiuDFLKTpYfokSzAeW1dlz67t+4/L2V9SryoRwZKOnZSn6xf7ziAOau8e6qJo2YvbnChZMamwPrD/M0AeU1dsxdewTHS2vwy1Y5bP6TFdwHfFO+OqlUY4ul+KJSHGFJL+04kx6p8bLm/OMm7eyXpdU2lTDPStROohdvMqB5slrQpniz3YPfYxP7tEScSb4XnvpxGxZsVqcXiDPq0dPNY0fye++aI5tuP7mOz5UUlNa4vre7eJzrt+kNUVmyMboEe7uRvGLMsHvC3ZPwISUYG/88zyqZ1JzfeFd/C5z1IPf+SGvLXwi5Q3mAkTQ5OOZxHnqf2gbI6oRnMl5Undpp9F8rkWzaUpBTvNlTY2+TLk+EJpoNqKix46Do9eFrKK/1Xf64Sz44b7PXbZWivdhdsB8+VYXvNgQp97ofzF+fjx82qgXkg99sxpRZK3HL/9a7fMnf+VP98pZC7osamXrZF1W1ssYu8deDo1zCMSVe25OotNqGDIX3k5bGLtEiNQ5/PzwaUwfyEajXSVkA256egFcnn4HbRnbAbSN5Ur8Nh+W8OV/fMgRvXdEHuZkJqKhVm/0kwT5MTICXYNKjv8L2L72wJDPdL1uPI5qJLsEOcC+SEEdxRRRpuVw4D7ldnVWy/Whg9Azu/aGk41her/WxIo9ixvnJ/fCEbRqG1ryOM2vegJ35P+ViF8fDkntjvNHTxh6vmGBLtBjwg0LDG/XKn35HD9qdTk1TTH2QBOSBokpVFakxr/6Jf31Vt307WNw7d5PmC0l6ARaW1+J0pafJY0IPHhy252RF0PpWaXWAMcCiMI0kWYx4RzSf2LykDi6tsqk0+wSz7/uqRWqcy9snOa7uezDBbMCDEzzzEw3ITcf5vVog3qRHtU2tOEj3WnvR5GVzCEiyGJEk9i1NfElZxPv4l23HcaAoCtMEi0SfYCcaht5TExrWORufOcbjKLJwDJlwuHmWzF+fj5X7tHPWSOkDJPdGl8auY2iXmYB4k96rvR3gaQL2+imUCkpqVBqgEn+DdKT+FpTWqIqASCMPu490CMHC3YzSpblnubiv1+WrUtcCPIpX0oIbo7G/t3Qf3l6y1+v2qlo74o16jxB7yWyh9NBRUlpt82lS0WJEJ66Q9GrVsHQer06SfeTjTXoPH/z3lvIRT67oTivdmzmip5ZWf08FcTQUbEiwN2GuGqSegHW4FQW5d+4mTP3gH81j7ZIpRnzoJRu7zSlg4V3D8cd9I1XBSolmzwenTENjX7rbM9HUsdJq5GZom4n+fGCkan394dOaI4G6ClP4ynMTLKrdQvannZnrsc8bv+/B67/vUbW9Mqm3y27d0KyZy3YX4sWfd+LlRd7LRVZaHa4XthLJNKP0klEiCfZ/n9dVJXB9MbFvS2x7ajzaZ/kfZd0pm+975+gOmNi3lavdqNdhw+ES1Nrl67vtGPcgap0Wj7vHdMTcm3nwUrMkLtjTFGalB8Z3BiCnPIhGSLA3YRhjWPnIaEzoziP/7PVIheoS7G7ujla7E3EmPZqnWFQau5a3x/rDnjnSH/tuq0dbjc3h4RMvkZ1sUWm6E9/5W1MLrVOw17ciUQCodAtKapUW52VPNRaDHinxRrRJj8eagz7yzHvB4RRwjSLbpuD2Qp/5227kPrwQlbV2TX91La8ViS9XH8bxshoY9DrcOLydSuD6gjFWp8nGnc+uH4iZk3vjvnGdVe2SR8wT329D7sML8ZeYlbJlahx0OoZ7xnZypUOQ7lul6Uiyw1d6GZFEAyTYmzg5KXEY2ZkPg+vjMSJpuJLwlh52p0JIGHXefdoBdYoACaWWJVFtdfg0uYzorJ5XWLLzpMc+kSjYq9wm+PzNpyO95Pq0ScX2Y/X3ZXfX8mvdfrs0Qth1vNzlEeONAc8txt1fyvlZpEne+kyON5SclDhc0sfzxfHkhd0BwFUb9+qP+EvsvnGdPPaV3CiVioM0IqkMVLH0MECCnXDZUO31CFG3u7k7armRKzV2LW2sVjRFKO3bzZIsqLE5sOcE9zsXBAHVPjR2AHhwfBeM7tLMta41oWirw9TiLtxCgbsd2OlnfVwpJ056gglHS6rrbQted0it5X/mZsOX2HWivE5beWF5Lb4XvXp+2lLgSs+cWk8beyDxNqI4t6dnNtKRnfl900eRlVS6V6tIYyeiGYMo2OujsUvh7FI0qJbPijJSVNBInCQJ04fny0E3W46WYuoH/2DszGWotTtgdTjhFGRvBS30OoZzejT3uh0Aat0E+4myGjys8EYJtY29vMaGc9/4S9VWV2TuI+d0QbvMBHQWTU/Se+D8N5f7OErru9UC64Wf5QIr7gFQ6Ynak9b/Ps8zW+Vts9e7lu8c3aFefQo0d43p6NGmdQ9d2rclNj0xDh2zZXOeNErxNjkcDZBgJ1wC2OZwusLB3e2u7qwVbbvdctxy9SgOU5oWlNXmJSQXxG/WqaNBJd/ksmq7S9j5EuwAcFk/77Zcp1PA+0vVSZyKKmpdQ3Wgbo0+0Lgn1QK4h89n18sl9a4bmqvafsWgNvjj/pEu84FOdPstqGdhZy1zl6T1z151WNWeFq+tefvKr59sMaBZsve8QKHgnrM7YuUjo10mPMlP3h3GmMeoRLK7a2W3jBZIsBMu75VHv92C3k//KpYJ862tFFXUYnC7dLQWg5C6t+CTUSMVJhGlxi7ZsF+fcoarrS7zR3mNzfVw+Uo6BfAHVInSZr72kOcEo/t3h9rGrqUNDmmf4XL7A+DhCeQ+TzF5gBzHoHSdfH/pPny/0XvQ1Yq9ni6sY2cuwyPzPdMVeDPFuJs7lN9f170TChhjyEmJw6PndYWOAY+e67surxKjXgeTQecxuR1NkGAnXAJY8rCosTkwZ/VhX4cg/3Q1WqbKkaXdWiRjy5PjvCbYenBCZ0wd2Brju8smE/cKPe6U1dhd+Uq03O7c6aVIMaA0a2iNPtwFeegFu9olc3jHTNfL6Y/7RuDnu4fj6sFt8fG1cpZD95dX5+ZJePx8nh5WcvGstTvwws87cfeX2kFXi7efwOIdJwAAC+4c5tL+iyutmv/zJC857t0LnJz2kjcm3Fw9uC32PX+u19/hjUSzAe8v3Y9J762se+cIhAQ74bKxS1jtTpWgcxeMDqeAE+U1aOmWhtfbwzOuWzayky14YWIvWIx6zJzcGzkpFmwvKMP//jkEo57h1pHtPYb3l737N6Z/zgsjx/sRiPTV9CF4+ByumZVUK+qkihrkJ9cOcKX9nTJL7Z8fasGutHN3zUnGg+NljbJdViK65iRDp2MY1bmZ1uEuJC1eGgFomXiUbM6XQ/B7tEzBh9f099hHcvcDtJNmAZ6mmPomdgsl7i9Ef5DMMasPFgc1H0+wIMFOQO8Wrn/5+ytVw1B3s0W1zQFB4GkC6mLrU+oMfgBwSZ9WOKN1KqqsDvz7u62wOQRYDHrMv+1M1X52p4D9hTysO95ct2CPM+kxKI8XLvlZkevjRDm3QbfPSnTV23THfXI12Ega9b1jO+Hnu4erEpopqUsouTw4xBHKukO+Baz7+YZ3zHRVqJK4YXgefr9vBG4anocLemmPwKRRnpT0bZki3e27V9azqHwEkqBw8wxmPp5gQYKd8LDd7i+sVPlH19rUQk/ynPBHy000GzT9s901PotRp8oM6U5d/tQSZ4hua68v3uMy4+w5UYEEkx6t0+O8TvqFw48dAG4a3q7OfQblpeO+sZ4+2ID8wquoteObdfl4aJ5sJ9fycnLPac8YwxUDZVt9+6wEjOyUhfZZiZhxXjePdAISbUX7/xtT+yA13ohliojhjtkhqB0cZApK5ZHPyXIS7EQUouVzrCw4sPWYOlWs5LfcGK8Bdy8Xi1EPgw/BXtfkqYSkkVbbHK5AG5uDR8NqeUDcKwrMcAl2i7HuR/Crm4fgTg33PQCuHDrHSqpx/9ebVNuOa6Qb0PKVVxalvn5Ynl+mi87Nk7D3uXMwoUdz9GqVyotNi2SH2SMmECjvxZIq/5LVRRIk2Amkekm/KnHlh6s02xuTS8NdsCdZDK4UwFr4ClDyxtJdhdhw+LQqn7t7ClnJne/OORs8jg8WHy8/4FpuiP1XSZfm3N10X6FnUNaRYs/C2DXi6GverXKh57HdsnHLCJ4Gtz65WiTh105RyWr2jYPqPVEZiXxwTT9Xrd1InRj2BQl2wme+bF80ZlLJ3STSNiMevlznE/w0xSjZebwcl7zzt6oCE2NM5bnTOt2zeHaweXrBdgBwCdPGYDLokGQ2aGqVxRrpfmtsDiRbDOjXVi6izhjD/eM64cvpgzG4XUa9+3CHIhipT5uGZWeMNPq1Tce3t/M5H/eKUNEACXai3kjVaW4f1fDoQnfBnhZvQk6KBSM6ZWHOTYM99vdnotYbVrcKTMqyaVlJZmQlmeudZjYQFFcGxnabmmDUzGiplT2z1u7QDPYy6HUNEuoAVJOv/qZRjgake0IrMV2kQ4KdaBAjO2chzUuOdH9wn1BNshhh0Ovw2fUDMaR9Boa0y8AVirTC/ibI0sLmcKoSkikrOyWYDBiYl44ML6HzwSRQaWF1jOFbRRWoa8X0v8c0IlKrrdqCPVA01rQUSSiDsPwtChMpkGAnVJzXK8flC+4ehCJRbXOoqiMFAve0vnOmD3b5nHdvkax1iN/U2p0wGmSBo8oLYtbDpNeFPKUAUL+ka744dEq2pQ/KS8dj53dDRoIJ290mvQFuYw+GVj2pfytXQe5Y5AG3ielIp+HjWyImGdkpy+Um1yk7EVuPcrfHKqvd5XLI0+g27taRBOnAvHS0So3T1CJT4ngZtjPb189EcO/YTnj1t92u9apah0rjH5Ar25fjTQYYdCxgQrY+1Cf/vb9M6NEceh1Dx+xEzfw81TaHX5449eWly/wrqBGt7NWYnI5kSGMnVEzs20olZO8Xc1grQ/R5Gt3G3TpWUZCO7JyFVyef4XW/c3vm1Om1487Vg9uq1lcfLFaZYpS+3HFGPYyG8GjstgC9TGZd3c+1LNmFE80GzXw0NTYHzDFkBw8VfVqn1b1TBEGCnVCh1zHXxGa80YAMcWJMSmtrtTtRXGlFY83DkiD1FZTUULSCapSmGADIE1309DommmJCr7FfMUg742B96ddWFjqSYE8wGzTd9GrszqDa2GONjY+PRbMks0dun0iHBDvhgTT/lWgxuASvFMDzh1id6Ms6koTVhSTYGzMp6g333De8Tf09X04fjL8fHu3aP5Qau8mgw80j2tWZB8Zf0hWT2FKJt7zMBBSU1uCkGKQkCAJ6PbkIm46UIC4IpphYJTXehDbp8ThZXqsZFxCp0H+Y8EByX5vQvblrAlUS7FIQ0X8u7dWo7wimYHcPm9f6nuxkC1qIASihNMXU2h2w2p1Iqmd9T18oPVEkjV2qAyvlVz9ZXutKp0sae/1Ishiw4XAJhr+0pM5iKJECCXbCg/656Vjx8GhMGtDaJRDnrD6C46U1rgjHPm0aZ3OU7ObpCYH3H1dq7JL3jsng3Q0vwaSHzSHUmUY4EFSKdU61asA2BsnTRZrgTk/gL2cprcIGReFw99EL4RtlJO3n/xyMiiLX5BVDaCKFU0v29o9XHMCKvUXYJdYidXdPrC93j+mI3Ix4VX72QKHU2DtmJ2HTkRKfwkyKvC2qqEWrtOBGolaIWrNWDdjG8NdDo7Bk50nXKGRgXrpq+3GFT/uPm4/hlUmx7cUSSJTBcc//tBO/bD2O+bcNDWOP6oZe3YRPlAWpJaEOAM0amIZAwmLUY/KANkEJaFGes7kYZSoVl9AiQ9RutULwA41Uli7Q5pDMRDMu799ac5vTKeCIIklXnltlJsI37krM+sMlXvaMHEhjJ3wiJZlSotexqIkwlOqC+spEKT247kWeg4HkTaQL4fU7XWXFR2Lisdk3DkLHGA4kCgbuZQCjAdLYCQC8ylFuhqcZIivJjJGds1RtWnm+I40nLuiGYR0y8cQF3evcNzFEgr3KandN0gZhztiD18T4gJ+2FLjahnbIDHuh6WjDpJF1NNLt7KSxEwCAWRol0iR8ZV2MVK4bmofrhub55e2SLE6OldcEz1dZEAR0e3yRK09NKDT25Dj+eD/2/TYAvMYpUX+kW+jO0R2Ql5mAe+duwomyGrSrR4rjUEMaO1Enk7zYbqMByavn7K7efcalCWL3EoCBQBAEfLriAG4Wa7ceFn2htVwyA01OSpxqXVlQg/CfGnFexGzQIU305hr9ylK89MvOcHbLJ40S7Iyxyxlj2xhjTsaYd5WPiGrO65UT7i40ivWPjfWou6pEKhhhD4Iv+ycrDuLJH7fj1+3qyVtvJecCSdecxiVPIzhSaUiLUa8q+PLOn/vC1aU6aazGvhXARADLAtAXgggK6QkmnxNgkudPoNLoKlHat5WEcvKUaBy1Co09WvLNN0qwC4KwQxCEXYHqDBG53O2l5mYsICUIq7U78eQP23DoVKVrW5XVjld+3eV6uOuLN/mtD5Fgf/6SngCAmZPJb72hSCY6s5vGHow8R4EicntGRBRXDmrT6LzokYqUJmFzfgk+/fsgpn28Gqv2nwIAvPnHXrz5x17MXZvfoHNLkabuhCr484pBbXDwxfNwSZ9WofnCGGREJ+4V1qtVCpSv4+EdM8PTIT+o0yuGMbYYgFZ44AxBEL7394sYY9MBTAeANm0Ck9WOCB3Nki1YeNdwzFuXH5ZqQ8FEmsiUNLODp6owedY/WPbAKJwWg5YaqmF7S1MQKo2daDwX9G6B0V2aIcFsgCAI+Pd5XfHZyoOujKeRSJ2CXRCEswPxRYIgzAIwCwD69+8fhQ50BABc2i/2ND/GGIx6hj93FaraD5yqRKUY2NTQilH92qZhf1GlR3u0BHgRHCkFBGMMNw5vh1+3nVC50s5ZfRhju2Wr6r+GEzLFEAS0i14UV9ai2soDUeIaKNirbQ60y5RD+KU88EI0BgcQLkwGHax2J5xOAYdPVeGR+Vtw++z14e6Wi8a6O17CGMsHMATAQsbYosB0iyDCz6kKqysVQUNzu1RZHYg36/Hm1D5476q+roRjJNajG5NBh/WHS9Du0Z9wqJiPyI6XeRYPDxeNijwVBOFbAN8GqC8EEVGcqrS6TDHOBmrYUq3YC3q3AMD92htzPiIySI2TU/le/dFqAJH1PyVTDEEAWPbAKNfy/ufPBQC8++c+V5qBhppOTlfaXMUvAIX/euTIAKIBaDkQMETOvAkJdoIA0EaRAE0ZFXqyrBaAnC/EH37YdAz3frURNocTh4or0TZdeW7+GQV51AgfaBVYN2okCwsXlASMIHxQIWbxq88w+645GwAA1w7NRY3NiZZpcs4WSWOPpGE7UX+0KmCZIii9Lwl2ghCZdXU/lyB3x+mnir3tWKlreX8hn1RrrpEmlwR7dKMp2EljJ4jIY5yPMn0OPwXx9mNlrmXJf715iizYB7fLwF97ilwl7IjoxGz0tGJ7UwrCAQl2gvADf23iRRVyeb3DYs4ZZWGLW0e0x3k9c5CbSeXpoplcjfKC+worYXM4Xamiw0n4e0AQUYC/pq1TBRkAAAspSURBVBgp3zoAnBAnXpU1M3U6RkI9BujRMgX/u2GQa71Pm1QAQFl18Iq11AcS7AShwc5nJqjW/bWJHymuQrZYQHvPSV78O8FEA+NYZFjHTAzITQMATBnAi9FsU5jiwgkJdoLQwD3S1N86r0UVtejZMlVctiLepA9JtSQiPHw4bQC+uWWIK1YhUtIKkCpBEH7wwDebkZFowugu2T73q7E5VAnDtLwniNghJc6I/rnpKK3iJphuEZLamjR2gvDCTcPzVOvXf7oWPZ9YhB0F3ofbtXYnzAYdpg7kqan91fSJ6CYl3ogeLZNdWSDDDQl2gqgH5bV2THznb6/ba2wOWIx6bDh8GgDPN0M0DSwGfYMrbQUaEuwE4QVv86Xtm3n3apE09vvHdQ5Sr4hIxaBnWLH3FPacKA93V0iwE4Q3Buala7brGMPvO054tAuCgBqbA2ajDt1bRoatlQgdO49zgf7499vC3BMS7AThlXHdm2P9Y2NxmVvVqM35pbjhs7UodjOz2J0CnAIfklsiKG8IERoMovdTs+TwV1EiwU4QPkhPMOGlS3tpbrO5pXyU6puajboGF+YgopfZNw4GgIgoj0eCnSDqQKdjLi8XJbU2tWCXimFbjHqYDfzRSk+IrcLfhHc6N09CZqLZVXUrnJBgJwg/kBI85Srytrt7QEiC3WzQQadjmDm5N767bWjoOkmEnfQEIwpKqwHwOZfHvtuKBZuPhbwfJNgJwg+kfC/jFRkgx85chms/We1a33eyAoActXpJn1aqAh5E7NOlebIrXXNJlQ2f/3MId3yxAZUhzvxIgp0g/OCOUR3wwPjOuH6YOmjpz12FruVrPuZCXjLDEE2P5DiDq5xijWJEV1YT2uRgdAcShB+0SI3D7aM6IDvZgreu6ONzX5o4bbokW4wor7FDEASs3HfK1V5ZG1q7e2TEvxJEFDHBrSDH9xuP4oBYVANQp+klmhZJFiPsTgFfr83Hg/M2u9qrrKE1xdAdSBD1xOBWSOHuLzeq1k160tibKslxXKRuOHJa1R5qjZ1MMQQRYLTKphFNgyQLT98rFVmRcI95CDZ0BxJEALmwdwt0bJYY7m4QYUIyw/2x86RqPdRZPkmwE0QAefyCbmCMCms0VZJFjR3gudrn3MSjUe0hFuxkYyeIRnDLiPY4VlKN9YdP45NrB0REODkRPpIVE+el1TZX9SyHM7SmGBLsBNEIHj6nS7i7QEQQ7q6uUmIwm4M0doKIeD69bgAOnaoKdzeICCMrST1ikzV2EuwEEfGM7Nws3F0gIhCLUY+8zAQcKKrE9UPzYNDxacxQ29hp8pQgCCKA1Irpm68Y1AYGfXhs7CTYCYIgAkjbDF46MS3e6LKxk8ZOEAQRxbx9ZV98cE1/ZCSaXTb2gpIafL7yYMj6QDZ2giCIAJKeYMLYbtkA4LKxv7VkLwBgaIdMtMsKfgAbaewEQRBBQq9XB6uFKmcMCXaCIIgg4Z6bvzJEWR5JsBMEQQQJo16HrU+Nd60Xltf62DtwkGAnCIIIIolmeSrzzjkbQvKdjRLsjLGXGWM7GWObGWPfMsZSA9UxgiAIomE0VmP/DUAPQRB6AdgN4JHGd4kgCCK2mHV1P9eyIATfp71Rgl0QhF8FQZBmA/4B0KrxXSIIgogtxnVvjn+f1xUA8NeeoqB/XyBt7NcD+DmA5yMIgogZLjqjJYZ3zETH7OD7sdcZoMQYWwygucamGYIgfC/uMwOAHcBsH+eZDmA6ALRp06ZBnSUIgohWspLM+PyGQSH5rjoFuyAIZ/vazhi7FsD5AMYIPoxHgiDMAjALAPr37x/axAkEQRBNiEalFGCMTQDwIIARgiBQcmqCIIgIoLE29rcAJAH4jTG2kTH2XgD6RBAEQTSCRmnsgiB0CFRHCIIgiMBAkacEQRAxBgl2giCIGIMEO0EQRIxBgp0gCCLGYKHIW+DxpYwVAjjUwMMzAQQ/JrfhRHr/gMjvI/WvcUR6/4DI72Ok9q+tIAhZde0UFsHeGBhjawVB6B/ufngj0vsHRH4fqX+NI9L7B0R+HyO9f3VBphiCIIgYgwQ7QRBEjBGNgn1WuDtQB5HePyDy+0j9axyR3j8g8vsY6f3zSdTZ2AmCIAjfRKPGThAEQfggqgQ7Y2wCY2wXY2wvY+zhMPWhNWNsCWNsO2NsG2PsbrE9nTH2G2Nsj/iZJrYzxtgbYp83M8b6hqifesbYBsbYAnE9jzG2SuzHV4wxk9huFtf3ittzQ9C3VMbYN2K93B2MsSGRdP0YY/eI/9utjLE5jDFLuK8fY+xjxthJxthWRVu9rxljbJq4/x7G2LQg989rTWTG2CNi/3YxxsYr2oPyjGv1T7HtPsaYwBjLFNdDfv0CjiAIUfEHQA9gH4B2AEwANgHoFoZ+5ADoKy4ngdd67QbgJQAPi+0PA/iPuHwueGUpBmAwgFUh6ue9AL4AsEBcnwtgirj8HoBbxeXbALwnLk8B8FUI+vYZgBvFZROA1Ei5fgBaAjgAIE5x3a4N9/UDcBaAvgC2Ktrqdc0ApAPYL36mictpQezfOAAGcfk/iv51E59fM4A88bnWB/MZ1+qf2N4awCLwuJrMcF2/gN8v4e5APf4xQwAsUqw/AuCRCOjX9wDGAtgFIEdsywGwS1x+H8BUxf6u/YLYp1YAfgcwGsAC8QYtUjxkrmsp3tRDxGWDuB8LYt9SRMHJ3Noj4vqBC/Yj4sNrEK/f+Ei4fgBy3QRnva4ZgKkA3le0q/YLdP/ctl0CYLa4rHp2pWsY7Gdcq38AvgHQG8BByII9LNcvkH/RZIqRHjiJfLEtbIjD7j4AVgHIFgShQNx0HEC2uByOfr8GXgDFKa5nACgR5MLjyj64+iduLxX3DxZ5AAoBfCKaij5kjCUgQq6fIAhHAfwXwGEABeDXYx0i5/opqe81C+czpKyJHBH9Y4xdBOCoIAib3DZFRP8aQzQJ9oiCMZYIYB6AfwmCUKbcJvDXeVjcjRhj5wM4KQjCunB8vx8YwIfE7wqC0AdAJbgZwUWYr18agIvAX0AtACQAmBCOvtSHcF6zumB+1EQONYyxeACPAng83H0JBtEk2I+C28MkWoltIYcxZgQX6rMFQZgvNp9gjOWI23MAnBTbQ93voQAuZIwdBPAluDnmdQCpjDGpsIqyD67+idtTAJwKYv/yAeQLgrBKXP8GXNBHyvU7G8ABQRAKBUGwAZgPfk0j5fopqe81C/kzxOSayFeKL59I6V978Jf3JvFZaQVgPWOseYT0r1FEk2BfA6Cj6J1gAp+o+iHUnWCMMQAfAdghCMKrik0/AJBmyaeB296l9mvEmfbBAEoVw+eAIwjCI4IgtBIEIRf8Gv0hCMKVAJYAuMxL/6R+XybuHzTNTxCE4wCOMMY6i01jAGxHhFw/cBPMYMZYvPi/lvoXEdfPjfpes0UAxjHG0sSRyTixLSgwuSbyhYK6JvIPAKaIHkV5ADoCWI0QPuOCIGwRBKGZIAi54rOSD+4UcRwRcv0aRbiN/PX5A5+t3g0+cz4jTH0YBj7k3Qxgo/h3Lrhd9XcAewAsBpAu7s8AvC32eQuA/iHs60jIXjHtwB+evQC+BmAW2y3i+l5xe7sQ9OsMAGvFa/gduIdBxFw/AE8B2AlgK4DPwb03wnr9AMwBt/nbwIXQDQ25ZuC27r3i33VB7t9ecJu09Jy8p9h/hti/XQDOUbQH5RnX6p/b9oOQJ09Dfv0C/UeRpwRBEDFGNJliCIIgCD8gwU4QBBFjkGAnCIKIMUiwEwRBxBgk2AmCIGIMEuwEQRAxBgl2giCIGIMEO0EQRIzx/xYZZlypQZzBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X_training_proj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimator 0: R2 = 0.994\n",
      "Estimator 1: R2 = 0.976\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "import numpy as np\n",
    "# Store estimators for later use\n",
    "estimators = []\n",
    "\n",
    "# Predictions for values. I'll just append here, and transpose at the end\n",
    "Y_pred = []\n",
    "Y_actual = []\n",
    "\n",
    "# Iterate over all principle components:\n",
    "for i in range(model.n_components_):\n",
    "\n",
    "    # Here, we lag X compared to Y. X will be one time period behind Y\n",
    "    X = X_training_proj[:-1, i]\n",
    "    Y = X_training_proj[1:, i]\n",
    "    Y_actual.append(Y)\n",
    "\n",
    "\n",
    "    X = sm.add_constant(X)\n",
    "\n",
    "    # Fit an OLS model and get a summary\n",
    "    est = sm.OLS(Y, X).fit()\n",
    "    estimators.append(est)\n",
    "    print(\"Estimator {}: R2 = {:.3f}\".format(i, est.rsquared))\n",
    "    Y_pred.append(est.predict(X))\n",
    "    \n",
    "Y_pred = np.array(Y_pred).transpose()\n",
    "Y_actual = np.array(Y_actual).transpose()  # Challenge - can you compute this using NumPy outside the loop?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5ef058f9b0>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXd4HNXV/793Zraqy5Jt2ZYt916RGy70YicvnQAhAQIEEhICJPmREJLgFPImJCF58yZ08kKIKQkl9BACNgHj3o2bXGRbtiyrWG1XW2bm/v6Y2TKaXe1K2q7zeR49nrlzZ+Z4pP3unXPPPYdxzkEQBEHkDkK6DSAIgiASCwk7QRBEjkHCThAEkWOQsBMEQeQYJOwEQRA5Bgk7QRBEjkHCThAEkWOQsBMEQeQYJOwEQRA5hpSOm5aVlfGqqqp03JogCCJr2bx5cxPnvDxWv7QIe1VVFTZt2pSOWxMEQWQtjLEj8fQjVwxBEESOQcJOEASRY5CwEwRB5Bgk7ARBEDkGCTtBEESOQcJOEASRY5CwEwRB5Bgk7DH4YPMenGrrTLcZBEEQcUPC3gMejwfnvbkAm/50c7pNIQiCiBsS9h7oaGsGAJztXZ1eQwiCIHoBCXsUOOdwd5wGAMhMy7ywat1GPPnD69Du9qTTNIIgiB5JS66YTOY/n67Bzh2bMeb4G6i1T8LXARTCBQAY+v4dOEfaj9dfeRyXfvmu9BpKEAQRBRL2biz913IsBQARgH+j4Vi+2gEAuPTgjwHEJ+xrdu7HiCFlGDW4NKF2EgRBRGPAu2I451BVHlffMrU5YrusqJHbZQWLXpmLI49c2Wf7CIIgesuAF/aXXvk7Xnvg8/D55R77eX0+OJjP1H7o+El8umIp/vmfNfD4lWA75xwnTxwFACzFlsQaTRAE0QMDSth37T+IF55/GpyHRuhf2HkbrhQ/gcfVhoNHj0c9193RFrG9efM/sFTciYs/XI6//ez6YPtzK5/FsSe/mDjjCYIg4qTfws4Yq2SMrWKM7WaMfcYYy9hZRcdfP4fr9n8bfr8/2CYwTeQV2Y/GZ74U9dyaA3sN+z6fdg0uh0bxNwjvorHNjcffWIUbDtyFheLuRJpPEAQRF4kYscsAvsM5nwJgAYBvMMamJOC6CWesUA8AUPxml4rs96FKPRb1XPWDBw37r7yyUmuXjdfa8r/X4fYtl/XXVIIgiD7Tb2HnnNdzzrfo2x0A9gAY3t/rJpq6upBo+2WzsPt9HqiMRT1/gW+tYf+6fXfh3Q9XQfEbY9ovklebzlV59OsSBEEkmoSGOzLGqgDMBrA+wrHbANwGACNHjkzkbeOiqylUKlDVXTGbdu5Ctd6m+H2o5Kd6dc1l/7kMbm4DYui2FxZwnwynlaJLCYJIPgmbPGWM5QN4BcDdnPP27sc5509wzqs559Xl5TGLbCccroZCEv2yFwCgbnwm2KZEGMUDQIt1WI/XdTJvzHs7mA+v/PprcVhJEATRfxIi7IwxCzRRX8k5fzUR14yXti4/2j3+mP3kMJeJImuhjcxZEjruiyzQXVe/YNhX7tpl2HfDbtj38sij8i/7/x7TRoIgiESQiKgYBuBpAHs45w/336TeseEXF2HtgxfH7Ofr6ghuv/7KcwAAZi8KtoWP2JtQjEPjbwIAlJaHpgtqr/oXxJJKNKM42NaGAsN9bMwYD3982IUAgIPWSTFtJAiCSASJGLEvAvBlAOcyxrbpP8sTcN24uEDcjIvETTH7KV2hnOq3n34YiqIYVpyePB3yHnXNvAljrv0t8O09cBSH3EZWpybiaphT3YborpgD0ngMv+3v2GuZAp/XjQ07KPyRIIjkk4iomE8454xzPoNzPkv/eScRxiUS2dth3Jd9UOWu4P5ZH2sx7Bsqb0XlZSsAUQIKNf96c6E22h46vAoAoIY9Njv3ol0oxsEZ3zHdM++87wIA/FIeJrOjmPfqwoT9fwiCIKIxYFaeqh5jFSRF9oNHiGfn1nygW9jjoG+vB1a0QbA5tT5hj83JvNg9/EoUL7rJdK2h868GANhdJ/prPkEQRNwMGGHnPrdhX5Zl+Go+MPVjVkfMa6nM+NiYxQmLzXweE0QAwHgheqoCgiCIRDNghB0+l2GXKzLOipCci1niEPZuj41ZnbBa7VF6EwRBpJacF3bOOY61uAG/W1tMpCPLoRDJ9YOvCW4LVmfMa5pG7FYnrN1G+m32EcHtlmvfCm5v/Gx//MYTBEH0gZwX9tUfvovKP1RAatoDDwsJO1cUuLgd64dcCy6F2qUILpXuCFwx7Iu2PAhSKH5967yHUfT9z4L7pZOWYPPk7wEA5v59bp//LwRBEPGQ88KeV/M6AGCeshVeZsOOCXcCAI7WN8AGH7hkB0RrsL9oy4t5TRs3hjg6CssM+5Kj0HxOeVVw2+V2m44TBEEkipwXdjhDoutjdngdgwEAU15fBompgMUB9HLE7rYPNewXDa0y7FudZmHPKwmlJjh1/FBcphMEQfSFnBf2Rp8luO0XbGCi5jKxM83HziQ7ykeFsgxLttg+9orbjVkTSrsJuy2vGN2pmrYAbiEfANB64mB8xhMEQfSBnBf2z9X9LrjtExymFLrM6sCYxVcH9y32/JjXtJYMQwdCXwD2fKOQ2/LNI3ZmscN9kxZe2dVYG5ftBEEQfSHnhT0cRbRD7RbPLlicYJINMtcehdUee8QOACdL54V2ui1ocuYXIRKlFaOhcIa2+oPoCEtc9vH6jdhzsDau+xIEQcQiZ4R9/ZatWL/FGJeuKqphXxbsphWogh6mGIiYscQp7ONuX4mWoYvhunWN6ZizwOyKAQDBYoPIOJY1P4sPf39LsH3Ju+ej4i9nxnVfgiCIWORM5Yf5b5ytbcwJFZ32eN0Il2lFcgI+o7CL+mTpydL5GNeyGnlOY7bGaDBbPkq/9nbEY5Y4Vq9e6nndsF/MXFF6EgRB9I6cEfZIeFwdBmFXJTuEsMRfQCi8cezXnofcsBeOwhL0lYbLXoJr978wpocSewRBEMkmZ1wxkfC4jRkdueTA1AtvMrQVFWtCzqx5kCrP6Nf9hsy6GGO+2HNK+qPzVwAA6qGlA/b7YxcJicb23Xux8ld3wOOTY3cmCGLAkNPC7u3qlh/G4kTB6GocdU4NtlVOXZJSm0YuuwcbBn8BeVyzbf3GdX2+lvTyDbi+ayVOHtqRKPMIgsgBclrYfV1Gf7rcqmVZLLv5xVCjkIZHkFeGQubGS7/5Jhb/6/PB5u6TvbEYodYBABjL6V8jQRC9JKsVQVUUU5sS1ubvJuyDvbUAAHtxRVLtioVQMAQAcE3nc4b2TpepBniPSHrOGr+3K0ZPgiAGElkt7LJs9i27XSG/uuzVhP3AOY8BAPLP+hYAQJAspvNSibVoSMT2ztNNvboOg1baL7xQN0EQRFZHxaiqWdg9rnYUFGpx5IpH82NLZWOAFW0YZuqdHpwlkd8YXO3NAMbFfR0BmutG8ZGwEwQRIstH7OaIEndnKI5d8WkuCmuENAEnyhbhyCxzndJUkF9kDKlshRY7v/2N/437GpzzYL4bEnaCIMLJ6hG7EsHHLntDKQNUv7YdKU3AsG+mr9720DEzDPunxXIUKx24yv9m3NdY//F7WKBvK+SKIQgijKwesfMII3Y5bPTKAyN2R+wc6yml2wKmNnvvnUTWTx4Kbqsk7ARBhJHVwi4rZh/7ng/+AlnWR/J+Tdhtjvjyv6QLpXR8cPtUy+m4zpHCwjSb9nyScJsIgsheslrYeQRXzGWul/HK079Ec1snoI9krXHkWE8nxZOWBrdr3ns8rnMUJga3L2p/OeE2EQSRvWS1sCuK5oo5lDfL0H5N/UP47E/XAXIXurgVLB2LkGKghD36qrnLgttWmy1SdxNWuRNHhMqE20UQRPaTeYrXCwILlJrHX4UjV7xlODbLtxlM9sDLrJFOTTvivQeB7+wHVrRBDM8GKdpjnss5xxDfMXQ4hifRQoIgspXsFnZ98pQJIpxF5YZjAlfB5C74kJnCDmcpUBBaqHRg3s8AAIq3M9oZQWTZjzLWhs6SKTH7EgQx8MhqYVf0yVMmWFA+cqLhGIOqj9hjj4AzgWFLvgQA4HEIuz8Q+WONXcaPIIiBR1YLO1c1VwwTRIAx1FR+IXgsj3kxr/ND+IUMHbF3w+HU66TGI+xer7YhWdEuaGX4IqVXIAhiYJLVwq7oYsZELUJk3E2PmPrILL7JyHTDRAlubsOEYy+h5tABwzGPX4E/LPOjT4/PFyQrdlfdCADweoy1XAmCGLhktbBzNeCK0RbQMtGc3Msvxi5Tlyk4mReDWAe8z10bbNuxbRPsD5bizf++Ptim+H0AACbawCyaq8lHwk4QhE5WC7uqGEfskfCLmR3DHolpvAZbNn0KAJDfvhcAcIUcSoEg+zRXjGCxweLQ/Ox1h/el2EqCIDKV3BB2IXrKG0XKnhH7EXFUcFtZ/WsACFZaAoAurx/rd+zGjg2rtH2fHxWTFwMAOo9uS6GlBEFkMlkt7IHJUyFM2PfP/oGxj5Q9I3YvQq4kJU8LhXSpof/bR4/dhfmvLsRFe38IAHC27EFZ1WRtew+tPiUIQiOrhT3kigmJ34RLv4dj9gmhTuUTup+WsaijFge3ueKD2+3CHHUXAKCJlWLuaW0Rlsi0AhtTvrACFj1dwkx5Bzy+vhfGJggid8g6Yeec4913/oEujy84eSoIRh97S/U9wW3L4PHIFiZd/5vgtuh3obO1ObjfJA2B2u3XZS82VmLau/695BpIEERWkHXCvuWTd7Fsw41Y8+z94BFG7AAghSX9sjqLUmpfvxAtqCuYqW36XVj795DQS74OWLk3uH9aKDGdrrpakm8jQRAZT0KEnTH2Z8bYKcbYrkRcrycE1ykAQHH9x8FcMUI3YZ9UfU5wu7gsvYWre8uIe1ajRpoAq+LCpaefDbaPY3UoYqGQxprixaZzvZ0k7ARBJG7E/gyAixN0rR4RLNqCo2rsQduJ/QDMI3bRUQT1e8fQeuXfUDF+dirMShyCAK+UD4vixkFxTNRuPCzap84xCQCgtp1IunkEQWQ+CRF2zvl/AKRkuBgQdgAYtP8lrS1CuKPgKETx9ItSYVLC8Yt5sKlu2FUXNhVdgC1O8+icW0LuphH3aOGPQmd9ymwkCCJzyTofO+eh7anCEQCA0MMCpWxElpwYrR6Fg7uhWvJh97ea+gj2wtCO1YlD4hg4uhpSaCVBEJlKyoSdMXYbY2wTY2xTY2Njn6+j+L2mtu4+9mxncNsOAEApOiBLeXAoHaY+Qr4xTbFbKoJNbk+JfQRBZDYpE3bO+ROc82rOeXV5eXnsE6JdR44k7Lk1Yi9BSKArTn0CAaqpj7XQ+AxlKQ9l/uM4cPhIsK2l3YXdB2qTZidBEJlJ1rliVNlnassrMIf+ZTMuhPznXWfcFrGPo2SoYV+2FKCctWPcszOwY9smvPyLG1D68DBM+etM+GjhEkEMKBIV7vgCgLUAJjLG6hhjtyTiupHwez2mtuKyYcm6XVo4PeLc4HbxpLMAhCYWGsWh8DErqiYZo314XmgEP+W1C3CV7/Xg/okjlCCMIAYSCXFOc86vS8R14mHhzh+Z2gQpt3zsU77yJ/h+9iqs8MPmyEdgvN3FnHDc/h6sg6tM50gVU4Hj+jYzum687U3JNZggiIwiq1wxPDwkRmd/4cI0WJJkRAtcN3+EfeNvxaCKUWD6iL3h6jeQH0HUAcBWEv2txd3ZlgwrCYLIULJK2E81mBfgjP9GbmY1LBk5FROv/y3AWLCtp/TEg7vVfA1H+nBFIk0jCCLDySphr6/Zampjttwv6MwCbyphIt+dssrowj6dHUq0SQRBZDBZJeyehhoAwN4RV6fZkuyhRRqCVp6PTpcrdmeCIHKCrBJ23qWtwKz8wkNptiS1CII2Upf68Ns6eeYKFLNOHNr2cYKtIggiU8kuYfe0QeUMzvwsSsWbAOxfehF7R1yN4WNn9Oq8RsdolI6dAwDoOkkhjwQxUMiqOEHB1wEXc6BAyK2VprEYNHoGBt36VK/OcX1zJ8rLRqK9VUvfoHoo3QBBDBSyStjLp52PuhPlmJxuQ7KAvLKRAID8fG1VLidhJ4gBQ1a5YsaedR0mX/fLdJuRseyf8FVTmyBJcHMbmNecSIwgiNwkq0bs4dQ7xkMW7ahMtyEZxIQv/gb7n2iEbcalGBXW3sXsEJSutNlFEERqyVphr/jepnSbkJFMuO1ZU5sfVjDZnGOHIIjcJKtcMUTf8DMLBNWc7pggiNyEhH0A4BNsEBUSdoIYKJCwDwBkZoVII3aCGDCQsA8AFMEGSabJU4IYKJCwDwDchaMxRd6N48fr0m0KQRApgIR9ACDnVQAApCcWpdkSgiBSAQn7AIBZ8wAAQ1hrmi0hCCIVkLAPAJgtL90mEASRQkjYBwCC1ZluEwiCSCEk7AMAJlrSbQJBECmEhH0AIIihzBEtp0+n0RKCIFIBCfsAgInW4Pbep80ZIAmCyC1I2AcAghQS9vKu2vQZQhBESiBhHwAUDx4e3LZwyvJIELkOCfsAYOTUhWhnBQAAK/el2RqCIJINCfsAwfvF1wAADUWz0mwJQRDJhoR9gFA+fi6aUQQvl7Dqzb9CVXm6TSIIIklkbQUlovd4mR0L2t4GNr+NTQ4rqs//QrpNIggiCdCIfQDhZfbgttpxMo2WEASRTEjYBxB+wRbcZpxcMQSRq5CwDyCK5cbgtqs9OStQ/7nyYbz/t0eScm2CIOKDhH0AMRghMT+79ncJueb2LetxOixNwcU1P8EFu+9LyLUJgugbJOxEn1EVFTPfuBC1f7o03aYQBBEGCTvRZxTFDwCYLW9PsyVELuH1efHh72/BkcM1AIDDB/dh354dabYquyBhH8CoitKv8xXZnyBLCCJEzaYPcW7ryxj1bDUAYPRz8zDxpSVptiq7IGEfQHTCWHDD5+tf3hg5hrC/8+i9ePex7/XrHpHgnKPdTTlvchVJiry8RlVUePwK3D45xRZlHyTsAwj3Te8b9r1dXf26niIbP2A8LISyvr4Oyxsex7KTj/XrHpH44NG7wX81Gl1uV8KvTSSe2tqDOHHyBDb+ZAk++u/LYvaXfd6I7Y0Nx3DwF/Phf7Ay0SbmHAlZecoYuxjA/wAQATzFOf9lIq5LJJbBVdMM+6dbGlFUWtbn663/xyO4MGx/1bMrcK6+XfH41D5fNxIfPPszVB1ciWLWifNZB8CAprZmOJxUzzXTKf2/xShkbgwDAC+gqByyqsImiRH7+7vag9tcVcH07cba3ZjGaxBsIKLS7xE7Y0wE8CcAywBMAXAdY2xKf69LJIdDy1/AzsKzAAB5z13Ur2tdePRhw/6Cw4/263o9cd7h32CsUI9BrCPY1ul2D8icN7Kiwier6TYjbgqZ27Av/rQYtp+XRu3fseOt4HZ7W2tw21W/P/HG5SiJcMXMA3CAc36Ic+4D8CIAin/LUMbMWw7fhP8CAJSztoRdd/uGVbAgtb7PqufmY/WfB17M/KZfLUfbz0al24x+w1XzlxPnHEs73g7ub3kzNFhQmg8Gtxsa6pNrXJaTCGEfDuBY2H6d3kZkKILFFrtTL5n5zmWwsP5F2fSFxceeSPk9080C31qUs/bYHRPI3j274PUlNpe/291pamtvM66IPufQQ8FtW1ttcPvwB08l1JZcI2WTp4yx2xhjmxhjmxobG2OfQCQN0epI6f12bPoEpxpOJOXaErLHJZGtbF3zHia9tAjrn7yr1+e266kr2mGeC6k/UoNtaz/A9vUfBttaThwGAGwou8LUv8hTF9pRKTKmJxIh7McBhE9Tj9DbDHDOn+CcV3POq8vLyxNwW6KvSPb8pFxX5ZFntWa89TmceuLKpNxTYAPPxx6grbU1dqd+8p8Xfo3Z72vpncc0fRijt5mThz4DABy3jjYdO31gA2a9dwVmvnt5sK39VC0AQBw23dC3gZWhVDkV3Ff9kSNn+opfVtDu7l+UWCaRCGHfCGA8Y2w0Y8wK4FoAbyTgukSSGDN9YVKuq/Tw5zRWPtDn67a0NEc95uMDt6TAzr/+v6TfY+m+nwe3S9XeJ45rO74XANBRbI6nGL7l16Y2f0cTAGDwtHOxvezzwfYmy3CUIjRxzjy9/1LzyzJWPf8QXC5zmOzGh69E60OzoCq58QbYb2HnnMsAvgngPQB7APyNc/5Zf69LJA97fjF2WrUSeR+99PuEXdfCFOyXJkQ85mF99+sfeOomAMBuxxzTsRPisD5fN9M4eaoBXV3xL7yyRlnIkyg8HuMI9jQrRltbfBPugWglf/MRAAArrQoeq5e0KbhhLPSFvfpvf0BHRzsUnxZBY3MWQpj0ueDxrjLjCF70Gu04WV+Hjx9choaGeihRIqU2vPN/OGf/g9j61x+Yjp3pXoWROIm1b/9fXP+/TCchPnbO+Tuc8wmc87Gc8wcTcU0iuaiS5mc/a88D6HSZJ7H6fF0WWWwURI5ZjocSn+afl7n5z9UrpHa+IJkMfWQCDv16acx+Pq49S8VakFR7Wk4eNewPRwNcD8/Gf35/I2p2bYx63vHa/aj/yTis/9dLEDpOoAMOTF3+dXxWsAjNX92Cih/uNp1z9u4fYcez34a3RfPi2h35mLzkcmwvvwQnLn8FtgnnBPs2owiSzyjsh9/6NZb4P8WQRydhy0PLItq1aMt3tWt31mH9qjfg1xfYNTc1hPX5dk+PJGuglacDFL8Umsza+fgt/brWYesEdMEKAFCYJdi+oSQ04ipDK063NPXp+l0WLebZ6Q+5Ak5II/CZbTZsqhurnrgX7R2JC91MBwEXwFR1X4/9fF4vrIHoIyWxUSrd8bo7TG3DWDOWtv4D3tfvjnpeU+1ODGdNGLT+l7C6GtAilMFZWIqp33kHg4aPjXqe3dOIpcefBABYnXmQbA7M/MZzGDbzfNgKBgX7NUtDMce9Brs2fRQ6OWzV81zP2h7/X9PaP8L8j76M1p9pfv99L95vOJ4L7hgS9gEKDxtZL2z/Z58Sgimc4Zg4EqN/sBENwlAAgCqErssl42h6z0s/jnqthpMnUHe0NuIxj7MCANBWNhsuONDGClH2nQ3wWItRpR7DOScex87nf9Rr+zOJzs74fMbujlA/lnRh10IqO7rlGAIALkZ3rfk7WwAA45RDmOb6FB3W+IIlzuhcHdy22Yx/O46CkuC226Zdb9pblwTb8k9tMvSX/dGfjZ1pOY7KWTs+/P0tEOyFhuMf/znx+Y1SDQn7AIULRtdIV4TRWSz8kHC8XMu6F3CJhI/YVYtREMSCwVGvNeSxyRjx55kRjykWLYpn1q2PIO+BehT96AisjjyAhf58JTl23phMXqVas+Wj2J0AuMO+APKbtibLHACAX/+bOL7sGVPEk1+MnspBcbUEt61MgWvwGb2+NxOM0pRXFBqx+5xDTP2ny8ZpvQO7NxttUpSIUVvntr4MzjkUzrDPqqXBmHr8b722N9MgYR+gDDn364b9rgiLRWIhQQFETcj9oibsPGzEzrr5gFV/3zIyMtWPDu6AaLEBjAH6h95fFAqh46znP+UTRw9C+Gkx1v9zZZ9sSDajP7ozZp9jdUfhbg+J5jTZ7KtOJP4uTdgtjkJsH/c1wzGLHP3vRXUbo2csFeaImF3FWlYhj+7CC8f73cOmtvyiUE4jVQidwzlHRwQ3XNmr1xj2O9uao4bGWloPoJM50TVRC7uUUryCOhmQsA9QRk1fgn3SxOD+rtd7Fx2jKgokpoKL2odMFrXRuSqERuwsLF7ewy2ArwMtTaewa/MnvboXU2XIESZly2eGT5L1nBnq1P4NAAD7tsyMemiStOie9ghuDwBoajiOyqemw/3mvSmzSfFowm7PK8DsL/8SXd85Ejw23bsFraejzJl0C0V0lFSYukz71itwf2s3jnzueUP7uiHXwpZvziNjtYe5ZnjIbbh/2xpsfeqbpv7hhdsBoEO3tRElpr7DXbvRyQow64rvYnvJhcjjXVD6Wasg3ZCwD2DCI1jOPvEkfN74R9R+fYEI04VcljRBCvfdW4tDmSVczAnB14mmRy7CtDc/17sJKsUfMapmbPX5we0Fp17Cmjf+HPUSTNCEn2XgSlWPx4MJsjZpWgg3PnnrWTTU12HjB68E++x5838AADP8xkpCyZzoU7zaqNyWp/mgHQXFUMK+QPe+H/lLUvS2oRlFwf2CshHmToIAZ+lwFA8eaWy3F/do02f5CyGIocGD690fY2mbednM8YIZhn13uxZaedIxLti2c/EfAQAVaIJLLAQYg6/yTFiYgs3vPtOjHZkOCfsAprv7QlXifwX1ByanJO1DpkrmEXvhsFBMu4c5IPo6MUE9BADoaAu5FGIhqH7IsMTst2jLPdHt7dJ88AI3jsSOHzmAde/+NW5bEs3at/+Cbf98xtC2eNO30PHk5zD345vh0e1eUve4oc+20osBAM2n6pAsuC7szvzQ5GLXt/aghWmjXoszsghbfG1oF0Ij45LBEYRdZ8ioiTj1pVXYWKItRmKOoqh98aMmTP3225h6zU+xdbCWcmCOL+RL33fBM9g+52c4KI6GKBszSno6NGF3j9DCSTu5A1XVoTc+j6Tdd/qFXwEA+Pf9K7odWQAJ+wDGJxl94HIvhF0JFEPQhVy16JNpYaMpmzN0/S4hD3PaPwjun/7DErS1xifujMuQWd/j4AFA1Re+FPkaUH8iFJ+d/+elWLD+G/26dn9YuPFOLNhmjsIYp9YCAFoa6tAaFmcdwDLnegBAw8FtyTPOp32p2B2h32N+aQXcX/wHAECNUkHL5m9HV9jfVl6h2f0RzuBxc8Bk7W1RjPJloR20AIKIvOIyzL7D/LYwcdHlmHnJt+AX7LAoRmH36cLuHK75+63wo6BoEFzcrh23ave152v/Lur4J44d2tOj3ZkMCfsARq4807Cv+I3Cvu7tZ7Bz3QeIhC/gipH0iSxLwBUjokkf0YmWkMj7RKdh8qoKJ1BfE19Uh91zCm6xh5FcHHBZs7eS16PiidAqxiKmiVdP4XHJonuI6XHB7Is+uOov2PHWH03t+WVaeiZPq1n042X1C7/F6pW/inqc+Trh4nYIovFLNa9Yi25SPJEbVUTXAAAgAElEQVQzTNqVDngthRGPRUOoWgQAGDJhXq/OC7Bt+g+D27Jgh0UxuhVtWzU3XWmlNq/EwAHG0CRqoZOKzfz35Wo8amrLFkjYBzBnXNN9YYZxBLZg412Y/k9zlj0AUHQhZIERulUfsXMFefdsReNNn0KUQtELPsmceMzbYR6xb13zT1Nbqb8BHc6RpnYAaLzFuAJy9/Z1aGszx4QHhD0aPm/qE0B1j0Q6Mdm8UGzJkT9iae0fDG2b5/0WJeXaZKu/tW9ZM7et+SfO3vdTnF3zi6h9mN+FLmY3tTsLtFHtgr3/HfG8fLUdsqUIm8Z83SC4PTH78rvhvnM3KidVx9UfAE7oaye2VFyDWVeG8uYoggUT5b3whpXYmy7vAgBUVE3B5mFfxNFLXwYAdFi0aBvuCE3Ybhx5q3adLJ5AJWEfwDBBwKmbPsVuhxZnrCg9F6cOR/FrfQPCzmyacAuKD47CEpRXTYUohUbscgRhdx1eb2or/7c5NayNe6BYIkeLlFdOwGGhKrg/5bWLsPOlFaZ+XO55RO7vxcRxonC7jGF69iHjovQ0IlrzUFg6GCdZOaz10Zf2R0ORZcx6PxQO2L3gxUfPPIDGByoh+t3wMHPKBps98u8CANa+9igGowWqaEX1Db80CG5PMEGAc1DvyjhIN72ObTN+hDm3G3Py++2aWG98/BvAiiJsW6v5y9uRB0EUccZtj2LsHC3cUtbDdJkz5C4qmqUVogkv0ZdtkLAPcAZXTUXnWG3iSoniM42E3M0Vwyy6APCQSEhhBT0Ci4zCEV1mN4KVmwXYzr3B3DaRaCo3vr5Xnvx3BIN7HrH7oxRQTiZel3FRWPgbTk9INifAGI7nT0eFu+cUBJFo7JYDxtvtS+2s2t+jnLWjyHUYXsE8Yg+nZsV0fPz7G4P7o3b8TrPRl3xRHDxyEmZd8V1T+/SvahPNi5v/DgCY9d7VAICaWd839VX0FbSBgQmgJSADANmTuBxKqYaEnQBELURRkeOfPFV1oQyGnkUII5fCsg+qERJWCYp2DTnsvi6xCB6PB4drtJWEXFVhhxewRBf2Obf+CbvsodWNo9Rjpj6xlt+nwxXTPReLqsb36m/P10aXsmMQhvFTaGk82av71u1Ybdjf9OrvIvYbrx6ER+w50dh4HMWSVm0ydeO/XsQwrn1Zy9b+zYn0h4Awd2fqhTebG6u0KJkhE+YHm+x6FJDq6f1q7EyBhJ0A01eL9ibc8cR2rehC0BUTPBKaIBUt4SNQc7y1qE9wrX/k1mDbaOUw9jy8DKNXngmvtws+nwci44AU/fVftFjhG3dxj/ZyxTgib6w/io/+eHtwX+7jqti+0tF+GmNfMRYTl73uKL01Nld9FZ1wYuRkzQ8dyHt1+pHzezjLjOIyrgxdXBMqP/fxn4wrTL3WyBEtG8aZXWZzPw09z7GXmVPjppJNw75karM7zW+Ncy6/C+137MSoqWHCnqd9KSk0YieyGaaP2GW/F21t2it0pELD4TiOrwEADB03W+sfYchuCXPFCF7z6EfQR9GLWl4ztM/2bQEAdLQ2wxOYYIxRzm92hFfyAKqiYmG9MVa9vbEOZzW9GNxXfKkV9lNH9pvaLFY7Plv8v8H9dYO/YDg+4bL7kL+iHtaAj1tfhzCWm99QeoLL5v9rix5SuaTxBUO73xZZ2HmUOQ9A82WXVUbOy58qzrj1D2j88upgIZZtxRdE7McEEYXdFkk59QVZ8JGwE1lMYPVo26vfRtHvKuH3edHe3nO2QYevBbut01ExRkucVDZGE3hpbChvthjmigmMzjefEQqvE9WexbSrvRldLu2LRrD1XM6ve9Ko8FW0ezeZfe7+bkIup9jH7nWZqxFNWXQJpp5/Q3BfLDEKTn63ePChS24KboeHTn6y8kFseOeZqPdmdeYJV88fF+GDR8yjcG6L7NYQJKPvPfz++bznN49UwAQR5WNnY+ekbwEAxt70WNznWqw2eLmFhJ3IbgIj9pk+La7c0+XCbn0JezRK/A1wOULVi8bMXIzOO3dj9iV3ROw/6upfYFvxBZh6XugVWVJ79nu7O1qDfmgxhrADwBFxVHC7K3xiMoLvuvsIPdUjdr/b+MW5yzY7+OV07NoPcPjyN3HG1fdh15mhHD7dv7yqpp+JdRO0qJN2PReK19uFxTUPYd6GyIWnt37wd8x3rQIAHLzk1WD7MDTivFPPmE+IssSfWYzC3toSqkeaSXVoz7juAfD7G1BQXBa7cxhuZsfC+r9i439fFLtzBkLCTkDolsLX7+0Cwl7Xu7tlFFlGOW+GUmAMT8uPEq62yzEX5SPGYdbdL8PuzMeWWT9BM4oxzFeLT1/+A3xcxNphN5jOG/z6F+F+6asAADGOAtyD7/4PNozWslZ2tjYG230uTUR3Ln4Um0Zp15v6r+sM56bax+53h6JGmtggOC4IxXtXTqrG6JlLIUgSpp5/Y6TTg4h2bXIzkM73ZAQXTziew+uC22PnnIedZz9t6rPfGsrGaCmJXHpQsBqFvXZr5IVsmUD3L6F46NLDPOd612Vl4Q0SdgKsW+3ME49dEVxODmijwHC63B3ahKYtdmm2zrtqMPHuNw1tcy67G0fyZqCIuXDmrh/ByhQwyY6Dl79l6FfCOjFZrQEASPbo+b8DOAqKUTxFm0g8siaUU9ujL+IZNGYWiudcFvFcJcFV72NRvUkbaW+s/ArKHjhkSGgWTvdRencCwh4oitGwa3XPN2bGuZBpSy+HC0bh8yy4G8eueR9byy7B7GVfjWKX8W/Gty/k7to277c925AFeMPi91uSmI8nWZCwE7A6jaFp09S9cJwOxUd7u4w+011Pa6NiHseCpvySwbDYzBOfqtgtZluyQ7JGn5CzRglh6874M7SFJ5MPPAV3wB3TsAcubkdF1SRYokzCqikW9gDTrvlpzD610hhsGnFTxGOinsfF62rHujeexLwdoSpVkcJXWbfEb0wQsWNEyD3WwMox/eyrUDl5HmZ/8zkIUQpmlwwbDwDYNv9htCMPI1pC5ehKxpiLjmcbpWooJXHrqexLLUDCTsBZaPY/zvJsCG4f3W1cIbqg7V1twxe7alE0VNE4SmRWBwRL9AU6Nkd8hZsDI9wS1oGtK7WUCUz1w8NsYIJgqMQDAOtG3qbZk2JXTIBIIXjdqfrhVlTfGnnOw6nnbelsPIIFW4yRQY31tab+HGb/t7MylDvnyJTbwYTYCdeqps2HfF89Zi27BUfsEzGCh2LpS4eO6uHMbCH0ZtPV1pxGO/oGCTuBvJKea1JOfd8cEwxoBTD6Cu8WVSE6iiD1IOzWOASwO+UNa7B36ydgii+Yz33QUGOkieDUcoRM3pC6uOs1L/46uB3L1RKLUVO0+Gv/SXMmwua6GvMJ+hfYvmWhUM8ZF9yAzUO01Zn5wybFfW/Jpr1huQvGBNt2n/MUCooHRTslazi57Cm06DnlvR19K8KeTkjYCZQNjZxgK0DUKAfe9yRJ3Ysh51eMDy80b8KR1/uVjBOUA5j0+uf0tL/6QipBwNaCs4N9AvlZnCx1rphFe38OAFhbYZ4w7i1Wmx0d3AHWZQ5P9baaV6QyuQvtcGLi/FAuciYImHnLn7DnvGcwZdF/9dqGcVc9ENyuOuPCXp+fiUycvwzq7R8DAGRX/LUDMgUSdqLXHBDHAgAmhX2ge023EXteyRAMHjYKNdYp2H3es6busXJ694Sg+g2l9bzOUHrcvEHD0Yp8dKDnBVDJQOxKzEiwQyiA5DXHxctus9gLigde2EztktWGyUsu79P9Bw0JDQwczvhcZtlAgf4mu2D3z9NsSe8hYSd6DYeA7fa5KC6PHAoX1zVEY0UkR0EpJKsd43+wFlOWXIYayyRsLlkePG6xmsUoXgTVDwVhRbZLq4LbNmchDuZX47RgrrOZbBhPTNFkFSKq298P7m8s0wRaaT1u6ivIXfCx+JKN9YX+upYyifAslvEWhckUcue3QCSErflLsG6MtlrPzyNPolm5B0oP2Rb7Qn6RUVjH378e427Q8pAfFqv6dW1B8UIJG7EXhE0W2vMKwAULxASJbG9gcSb9isUIXh/c3meZjDlfexKnUQh78y5TX0Hxws/6/iUZjS1FF6CeDU74dTOFA0/2322WSkjYCQPymAsg5Gkie0SqCra7O0O5w22qB4rYP2FneqjkPstk7ChYCrvDHKdeVFKOrfN+g9I73uvVtdePNMZeW5QuKCz0hjBxXsgPbM8rhCpYIKVB2AUef5rkeGkd83mIkgUnrFWw+c1L4kXFA1+MVLx9Yc49L6PigQiTtTnCYPfBdJvQK0jYCQNzLvkGBH2lHgfDulFatj9PWLUfG3rOjx4PgRj4lhHnYcZ33ozab/byr6Jo0NBeXXviJcawvyn+XYYRe3gBEIezAFy0QEI6RuyJuee2M0Ol8wJRPn4pDzbVHI4qqR7IQuJH7LlOfeHMdJvQK0jYCQOiJEHQF/H4BTuEgiHatp5Lxef1ogT9L6LA9DwxLM7iEr26tmheVKN2Wyl5nGlfFqIkAUJ6hF2svikh1xk5M5R4zZKvCbss5SFfaTP1taheEvZe0PbNvWhGESxyduVmJ2EnTDB92blfygOTNBGQdWHftVqrSlPd9Hr/7qEEaqYmXtilCKslVWacrLXe+i4abtayHHLBAksKXTFebsHaii9hxrnXxO4cByVhk9j2Al3Yi0ZjCJrRVK+tmuSqiuYVIzFB3g9FTLwrJlcpKqvAKctwOLzNOFHb+2pV6YKEnTDhKNZGs75xy4OrQWU9DW5gifmWWbGXwvcEU3X/chKEPdzVEkARjG3lw8dgyEgtZ3gqXTFejxs25o8rz068hEeiOPRVxLYR2gTx/ne0/O5NJ49iELQRPAl77/BJ+Zgk78GwZ+ahy5UdqXxJ2AkTk+ZdiJNf2YC5V94Npsebn1j9JE4dr0VX/V4AwOApS/t1D9WuxaVbCxK/SlEKE3Y31944uBA554lmRB6sTIGnq+8pEuLF3aGJK4uS57yveLn2f7bpRSJsRZoL7cxjWqHnup0fBft2f3sheia8EPu2134LV4fZxZVpkLATZhjD0FETAQCiXgVpYcML6Hz6Uiw8oNXHdBb0fcEQAMz68i+xcfqPMevCxIeRiWE+9jqLlrekJzETC7Q3lNMpyOLn7tAWDQn2xC7k6bxjGzbO+CmGVmoraacsMJYK9LacCG7PaI1Q7JuIihxWr3fhgYdx9A89l2HMBEjYiR4Rwka/Y9Ta4HbZ0Mp+XdfuyMPcK7+TlAUt4dd02bTY6snta6L2t+qj286WhoTb0p1AXVMxRqm/3jJoyAjMvSJycQ1VUYDW2uD+SbEiYj8iMt0LsU/2706TJfHTw/spQQAVE6qBD41tCmcQs2SFIdfHLnksevZGa55WJcjbaV6Wn2hUrhdtYLEzKCaK1uaTWNCgJf3add5fMHTsrJTdOyeQsm9OIjs+nUTS2eY8E3XMPJIrG1qJ7fa5hjYxg0qfRWPdxHuxyzYbI774h5h9bXqCMZ+7/2GcPeHubAsW9BDE5H/0Ns3R6svWrAoV8p625FKUDcuFtLopRDS78VwdPdcETjc0YicAALPufTfqMRb1SOay4Lr7AdwPfxxFqh2BEMEISbMSBVdVOH8zEk6muX1SMWK35GtvIvP3/AIAcOCytzAu6XfNQfTSkGtH3AypfDzmbr0PLSePIK8gcj3YTIBG7ERM/DOuT7cJfSaQPGyrc2HUPla9whNPQrENrqpY/8LPsf03nwMADOOaHz/SIqpEUzB4tGF/3KwlSb9nTiLrpSElG6wFWjhp5cqlWPfEnWk0qmf6JeyMsasZY58xxlTGWHWijCIyizOWfyXdJvSL09/Yiyl3vRb1uKjH6sdT6q+3rH/xQczf92vMcn9qaI+nSlF/GTNtftLvMRBgsvbWxywOSGFF1Rec+Eu6TIpJf0fsuwBcAeA/CbCFIJJCSXkFbBHqrgYICDvUxAt70eF3IrZ3rz1KZC5M0YVdssMSR1H1TKBff12c8z2c8+xZZ0v0mbWVt6XbhKRhCQi734sNf/wKjh8Mpbt1d7Zh/ZN3wetxRzk7FpFnKFhPC6YSyPppKwAAG8/4VUrul4sIurALFjusYVlIA4vCMhEaNhBxMW75ndgrTky3GUkh4Ie3NW7HvKZXgb9egT1rtIyTO57/IeYffwbb3vhjD1eIjlWN/IXAUhAVAwDzr7oHWNGGuf/1tZTcLxdhE7Q0z4MmLQTC3rT2ODPX+xxz2MAY+zeASHlT7+ecx50JijF2G4DbAGDkyJ5rbBKZR3nFSJT/aAM2vP4obEWDkV1JTHsmsFJV1Edmw3kDhr//JZyoWAuxS69Q38coFouexbI7QopG7ET/qV5+M1xLr8aY/AJwVcXa8f8PlQeeg8Aj/24zgZjDBs75+ZzzaRF+epXej3P+BOe8mnNeXV5e3neLibQy79KvY+bZV6bbjITCBAE+LmKGZ6OhvenoHgh6RIRo65tv9WRRlK9Alo1BpAOXvHxt9SkTBCy8/odol8ohhX1pb/z7Q2huOJou80yQK4YgAFiZuUydr/0UREUX9rBoiN4gKl04xkJpdY+wEdoGz/xFXkR0FEGCqPqhyn4cP7QHcz97EKeevi7dZgXpb7jj5YyxOgALAbzNGOtdDTOCyGDkjkZIAWHvY24XSXbDKziwad7vsfnMR+Cy6MnTAqkFiKxEFayYKO+F8PMynD72GQCgwN+UZqtC9MvRxzl/DUD0AGGCyEL8XISFKRA662FVtMlP3kchtihu+EQnqvW1AJ9tehRAWM4YIivxWYsAfd3StFW3AAAEZM5bGLliCALA8Rs3BLfFH2sjr3n1z8MeqBvKza6aeHDKbfBK4dkBtY8cU0nYsxnFbq4jkDmyTsJOEACA4aNDoZzhCboGKVpUDFfiF+Itbz6Krb+9DH6fB0OVengKQkm3uD5pqmaSChC9hutFw8ORM6iACcVcEUQPBNL98l6M2Ods/j4AoGbXOoxnfrDi8PBePRqmj28ARGYQqQIWCTtBZCBbFj8Bf1c7ImVY4Wp8Qnxw+8cYq2+31e0BAFhLhgePq/oCF06umKxGtJuFXSFhJ4jMY87510Q9Fq+wNx/aFhR2+VQNACCvPDRi7xy2GDi0GQUVYyOcTWQLgr5aOZx8JXNqoZKPnSDiIN6oGLU9VF5Paj8CACgaHCojuOD6FTh240aMnjQ7sQYSKaVwxFRT2zDeAL8v8amf+wIJO0FEoItbjQ1KfCN2oa02uO3wnAIAOAtDE22CKKBy9IR+20ekl/EzFmDHBS8F93daZgAAOk5nRiw7CTtBRIDdd8zYEOdkZ57rGBqghcIN89VqbflFiTSNyBBmLLoYO6xa/diuyVcBAE7sjl40PZWQsBNEBOx2YwFjNc7JTqf/NE44tdDJErTDxe0QxdQVriZSy6g7XsPuz78BKU9bUTxmdWZUVSJhJ4g4mLf9h9j5fuyKORbugyKG0g+4mDOZZhFppqi4FFOqz8LYeZ8HABy0TU6zRRok7AQRhbVDv2jYn77mTmBFEQ5v/yjqOVbuhSrasHbQ5do+Mje1K5E4ikpKUSOOgyqao2XSAQk7QUSBRcnAOPzVK6KeY4EfXLJjyOmtAIBidCbFNiLzkAUrRNWbbjMAkLATRA9EFvZaaXTUM2zcBy7a0Lboh8kyishQOBMxzbsNRz5bn25TSNgJIhqWCedHbJcgY9d7T5vauarCBh+4ZMewiXOSbR6RYVT6DgIAXK9/N82WkLATRFTOOO9qtNx9FKsHXWtoH6McxrS138bpxhOGdln2QWQckGywZUk1eyJxuKBNlHfaK9JsCQk7QfRIaXERlt7xaMRjst84Merp0gtXWxywOUjYBxrua18FAKiOkjRbQsJOEDERRAHrSy8xtfs8Xd32NWEXLHbY7NrorR5lyTeQyAjGTJqJJhSD+d3pNoWEnSDiQfJr0S21gZqlAGRPOxCWHMzv1T7QzOKAIIpYN+9/oX6FqkUOJNqEYjjdxwFoieN2/e4y7HrtoZTbQcJOEHHgs2qv18crQhOqo/52IQ7+8szgftOBzQC0ETsALFh+A4aPGpdCK4l0czpvHEq9dQCAtpZGTGtbhWnbH4SrvSWldpCwE0QcjL5iBVaP/z7GX3qvoX2sb29we9p/vgYAEPpY+JrIfmRbAZxce3PzekJrGDpJ2Aki8xg6fCTOvv4+DB4yHOsWPdVjX8lKaQQGKtxahALuAldVnFz7t2C7p7MjpXaQsBNEL6k++zLD/o5XfoUdj90c3LfmUTbHAYu9EBJTsfWV32DmZ78KNvu62lNqBlVQIoheIlmMJdBm7PyFYV+0GjNDEgMHwVEMAJDrtxvafW4asRNEViPZC9JtApEmRKf2tpavR8YEUP2pzSFDwk4QCeSj4itQNXFmus0g0oRFz8s+xaMlgauBVhaRq3JK7SBhJ4gEMu36X4Axlm4ziDRhyysObrciH+zyxwAAqpJaYScfO0H0gw9H3gl/52kMbd2Mgi+vxJjy9OcJIdKHvSCUTqAYnWgRtfkYrvhTagcJO0H0g3Nv/nm6TSAyCGu35G+CpAk7jdgJIgvYev4L6GioxdJ0G0JkFKVDRxn2BVGTWE7CThCZz+zFy9NtApGB2O121LEKjOD1+KjiFowTrQBo8pQgCCKrsXAtnXPVOTdCtKRnxE7CThAEkUCO2icAAApKhkAKTJ7SiJ0gCCJ7GXv789hw7osoLR8KQdJG7GLLQWx/7l4gSoH0REM+doIgiARSWlKKeUuXAQBESfOxV9c/DwA4tv86VE6cnXQbaMROEASRJETJOHb2pihnDAk7QRBEkuhe+zZVWR5J2AmCIJKExWJF5/ebgvu+ptqU3JeEnSAIIonk20NpnmdtuT8l9+yXsDPGfs0Y28sY28EYe40xVhz7LIIgCCKZ9HfE/j6AaZzzGQD2A7iv/yYRBEHkFhvPfi64zVU16ffrl7Bzzv/FOQ9E3q8DMKL/JhEEQeQWc8++BGvHfxcAsPazg0m/XyLj2G8G8FICr0cQBJEzTFz4efy5yYflQ0tid+4njMdYCcUY+zeAoREO3c85f13vcz+AagBX8CgXZIzdBuA2ABg5cuQZR44c6Y/dBEEQAw7G2GbOeXWsfjFH7Jzz82Pc6CYAnwdwXjRR16/zBIAnAKC6ujo162oJgiAGIP1yxTDGLgZwL4CzOOfuxJhEEARB9If+RsX8EUABgPcZY9sYY48lwCaCIAiiH/RrxM45H5coQwiCIIjEQCtPCYIgcgwSdoIgiByDhJ0gCCLHIGEnCILIMWIuUErKTRlrBNDXFUplAJpi9kofmW4fkPk2kn39I9PtAzLfxky1bxTnvDxWp7QIe39gjG2KZ+VVush0+4DMt5Hs6x+Zbh+Q+TZmun2xIFcMQRBEjkHCThAEkWNko7A/kW4DYpDp9gGZbyPZ1z8y3T4g823MdPt6JOt87ARBEETPZOOInSAIguiBrBJ2xtjFjLF9jLEDjLHvp8mGSsbYKsbYbsbYZ4yxu/T2UsbY+4yxGv3fEr2dMcb+oNu8gzE2J0V2ioyxrYyxt/T90Yyx9bodLzHGrHq7Td8/oB+vSoFtxYyxl/V6uXsYYwsz6fkxxu7Rf7e7GGMvMMbs6X5+jLE/M8ZOMcZ2hbX1+pkxxm7U+9cwxm5Msn1RayIzxu7T7dvHGLsorD0pn/FI9oUd+w5jjDPGyvT9lD+/hMM5z4ofACKAgwDGALAC2A5gShrsqAAwR98ugFbrdQqAhwB8X2//PoBf6dvLAbwLgAFYAGB9iuz8NoDnAbyl7/8NwLX69mMAvq5v3wHgMX37WgAvpcC2ZwHcqm9bARRnyvMDMBzAYQCOsOd2U7qfH4ClAOYA2BXW1qtnBqAUwCH93xJ9uySJ9l0IQNK3fxVm3xT982sDMFr/XIvJ/IxHsk9vrwTwHrR1NWXpen4J/3tJtwG9+MUsBPBe2P59AO7LALteB3ABgH0AKvS2CgD79O3HAVwX1j/YL4k2jQDwAYBzAbyl/4E2hX3Igs9S/6NeqG9Lej+WRNuKdOFk3doz4vlBE/Zj+odX0p/fRZnw/ABUdRPOXj0zANcBeDys3dAv0fZ1O3Y5gJX6tuGzG3iGyf6MR7IPwMsAZgKoRUjY0/L8EvmTTa6YwAcuQJ3eljb01+7ZANYDGMI5r9cPnQQwRN9Oh92/h1YAJVAOfRCAVh4qPB5uQ9A+/Xib3j9ZjAbQCOD/dFfRU4yxPGTI8+OcHwfwGwBHAdRDex6bkTnPL5zePrN0foZuhjYKRg92pNQ+xtilAI5zzrd3O5QR9vWHbBL2jIIxlg/gFQB3c87bw49x7es8LeFGjLHPAzjFOd+cjvvHgQTtlfhRzvlsAC5oboQgaX5+JQAuhfYFNAxAHoCL02FLb0jnM4sF02oiywBWptuWAIwxJ4AfAPhxum1JBtkk7Meh+cMCjNDbUg5jzAJN1Fdyzl/VmxsYYxX68QoAp/T2VNu9CMAljLFaAC9Cc8f8D4BixligsEq4DUH79ONFAJqTaF8dgDrO+Xp9/2VoQp8pz+98AIc5542ccz+AV6E900x5fuH09pml/DPEQjWRr9e/fDLFvrHQvry365+VEQC2MMaGZoh9/SKbhH0jgPF6dIIV2kTVG6k2gjHGADwNYA/n/OGwQ28ACMyS3wjN9x5ov0GfaV8AoC3s9TnhcM7v45yP4JxXQXtGH3LOrwewCsBVUewL2H2V3j9pIz/O+UkAxxhjE/Wm8wDsRoY8P2gumAWMMaf+uw7YlxHPrxu9fWbvAbiQMVaiv5lcqLclBRaqiXwJN9ZEfgPAtXpE0WgA4wFsQAo/45zznZzzwZzzKv2zUgctKOIkMuT59Yt0O/l78wNttno/tJnz+9Nkw2Jor7w7AGzTf5ZD86t+AKAGwL8BlOr9GYA/6TbvBFCdQlvPRigqZhketEoAAACySURBVAy0D88BAH8HYNPb7fr+Af34mBTYNQvAJv0Z/gNahEHGPD8APwGwF8AuAM9Bi95I6/MD8AI0n78fmgjd0pdnBs3XfUD/+UqS7TsAzScd+Jw8Ftb/ft2+fQCWhbUn5TMeyb5ux2sRmjxN+fNL9A+tPCUIgsgxsskVQxAEQcQBCTtBEESOQcJOEASRY5CwEwRB5Bgk7ARBEDkGCTtBEESOQcJOEASRY5CwEwRB5Bj/HzbDt81FIE9HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(X_training_proj[:,0])\n",
    "plt.plot(Y_pred[:,0])  # Very close!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.82854666, -0.21165555],\n",
       "       [ 0.85007914, -0.21861267],\n",
       "       [ 0.87956395, -0.19715205],\n",
       "       [ 0.9049668 , -0.19738715],\n",
       "       [ 0.83715122, -0.19311534]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y_actual[:5]  # Note: this is projected onto the components of the PCA!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_actual_original_space = model.inverse_transform(Y_actual)\n",
    "Y_pred_original_space = model.inverse_transform(Y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.5026558 , 2.58232715, 2.81036324, 3.24342046],\n",
       "       [2.51533851, 2.59451078, 2.82154773, 3.25223322],\n",
       "       [2.51207196, 2.59955972, 2.83794003, 3.2842495 ],\n",
       "       [2.52172571, 2.61093996, 2.85141202, 3.29977885],\n",
       "       [2.49352747, 2.57919097, 2.81557381, 3.26066766]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y_actual_original_space[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.020585769074154"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute sum of squared error:\n",
    "np.sum((Y_pred_original_space - Y_actual_original_space)** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's see how it compares to the testing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.07706719,  0.07938362]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testing_proj = model.transform(testing)\n",
    "testing_proj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.07938362])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testing_proj[:,i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.        , -0.19375164])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-2.050442147291337, 0.07853803293296605]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#\n",
    "testing_prediction = []\n",
    "for i in range(model.n_components_):\n",
    "    # Create a data row - remember our estimators have a constant, so we need that\n",
    "    row = [1, X_training_proj[-1, i]] # this is the last row of the train data (used to predict the test data)\n",
    "    # Reshape, so it is a (1, 2), rather than just a (2,) array. This is the shape statsmodels requires\n",
    "    row = np.reshape(row, (1, 2))\n",
    "    # Predict this row\n",
    "    p = est.predict(row)\n",
    "    # Statsmodels produces a prediction for each row - we just want the value.\n",
    "    # Transofrm the (1, 1) result into just a (1,), and append to our predictions\n",
    "    \n",
    "    # Potential error here:\n",
    "    testing_prediction.append(p[0])  # If this errors, try p[0][0]\n",
    "\n",
    "# testing_prediction now contains the predicted, projected values for the last row\n",
    "testing_prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.23307661, 1.19361555, 1.29270013, 1.65309461])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Unproject those value, i.e. transform back to the original space\n",
    "predictions = model.inverse_transform(testing_prediction)\n",
    "predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>2Year</th>\n",
       "      <th>3Year</th>\n",
       "      <th>5Year</th>\n",
       "      <th>10Year</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2019-05-16</th>\n",
       "      <td>1.22</td>\n",
       "      <td>1.19</td>\n",
       "      <td>1.27</td>\n",
       "      <td>1.64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            2Year  3Year  5Year  10Year\n",
       "Date                                   \n",
       "2019-05-16   1.22   1.19   1.27    1.64"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testing  # Actual values for the last row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0008708345194396898"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sum of squared error\n",
    "np.sum((predictions - testing.values)** 2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
